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ABSTRACT 


The article presents a detailed discussion of the theoretical models of four different fluid dynamic devices: 
an ideal propulsor, an ideal fluid brake, an ideal screw propeller and an ideal turbine. The four models 
are presented with all relevant mathematical formulae regarding the forces, the power and the efficiency. 
It is demonstrated that the application of the model of an ideal optimum fluid brake according to the Betz 
theorem for determination of the maximum effectiveness coefficient of an axial wind turbine is not correct. 
In the case of a turbine the inclusion of important rotational flow losses may increase the maximum value 
of the turbine effectiveness coefficient above the level defined by Betz. Therefore the model of an ideal 
turbine should be an inversion of the model of an ideal screw propeller. This conclusion is supported by 
numerical calculations. It may influence the design procedures of wind turbines and may lead to increase 
in their efficiency. 


Key words: ideal propeller, ideal wind turbine, Betz theorem, maximum efficiency 


INTRODUCTION 


A review of literature concerning the theory and design 
of wind turbines shows serious deficiencies in proper 
understanding and correct interpretation of the results of 
calculations based on modeling the flow of ideal fluid in an 
unbounded domain by means of principles of conservation of 
momentum, moment of momentum and energy. An example 
of this may be the maximum value of the efficiency of an ideal 
wind turbine, determined by Betz [1] and quoted as an absolute 
upper limit even in the most recent publications. 

Starting from the 1960s the models of an ideal propulsor 
and an ideal propeller were used at Institute of Fluid Flow 
Machinery, e.g. for determination of the tunnel wall influence 
on the hydrodynamic characteristics of propellers [3] and for 
determination ofthe hydrodynamic characteristics of propellers 
by means of the vortex theory. The experience acquired over 
many years enables a thorough analysis of the applicability of 
the above mentioned ideal models to problems of operation of 
real fluid flow machines. 

The principles of conservation of momentum, moment 
of momentum and energy in connection with the appropriate 
vortex models enable the determination of maximum achievable 
performance of certain fluid flow machines, such as: 

- an ideal propulsor, 

- an ideal screw propeller, 

- an ideal fluid brake, 

- an ideal axial wind turbine. 


In the cases of an ideal propulsor and ideal fluid brake the 
required results may be obtained without resorting to vortex 
models. However, if the information about flow details in the 
close vicinity of the disc is required, the application of the 
vortex model (or the equivalent accelerating or decelerating 
actuator disc) is necessary. In the cases of an ideal screw 
propeller or an ideal turbine the application of vortex models 
is necessary at every stage of calculations. The detailed 
presentation and discussion of the computational models for 
all four above mentioned machines is given below. 


AN IDEAL PROPULSOR 


An ideal propulsor is the simplest form of the fluid 
propulsor, in which a surface perpendicular to the flow (not 
necessarily of a circular shape) covered with pressure dipoles 
accelerates the flow both in front and behind the propulsor (cf. 
Fig. 1). At the propulsor surface itself a jump in fluid pressure 
is created. Consequently, a propulsor stream is generated, in 
which at a certain distance (typically about 3 diameters of the 
propulsor disc) the velocity achieves constant values at the 
cross-sections of the stream both behind (section 3 in Fig. 1), 
in front of the propulsor (section 1) and on the cylindrical side 
surface of the stream. The sections 1 and 3 are referred to as 
"far in front" and "far behind" the propulsor. In reality the 
closest to the ideal propulsor 15 the magneto-hydrodynamic 
propulsor, used as the "silent" propulsor on some submarines. 
In these propulsors a strong electromagnetic field generates 
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Fig. 1. The stream of an ideal propulsor 


motion of particles of salt water. Strictly ideal case of such 
a propulsor would require elimination of the boundary layer 
on the propulsor channel walls. Within the model of an ideal 
propulsor two cases may be distinguished: 

- anoptimum ideal propulsor, 

- anon-optimum ideal propulsor. 


The above cases depend strictly upon the uniformity of 
the velocity field inside the propulsor stream. Namely, every 
deviation from the uniform velocity field leads to the reduction 
of efficiency. Further discussion is limited to the optimum 
case only. 

In an optimum ideal propeller the axial force acting on the 
fluid is equal to the multiple of pressure jump and the area of 
the propulsor surface: 


2,5 " 
T-nR5(p, -p;) (1) 

In the case of the optimum ideal propulsor the distribution 

of the pressure dipoles over the propulsor surface is uniform 
and the pressure jump Ap may be univocally determined from 


the Bernoulli equation for any streamline in front and behind 
the propulsor (Fig. 1): 


Иру? +P =4P(Vo t+ V) +р; 08а) 


2 2 п 
Ир, + Va) + Po =K4P(Vot+Vx2) +p, (2b) 
After subtracting (2b) from (2a) the pressure jump at the 
propulsor surface Ap is obtained: 


" X 1 
Ap =р,-р, -pVom, ms) (3) 


Substitution to formula (1) yields the thrust of the 
propulsor: 


4 
3 


T= onim, eom. (4) 


where: 
па = V,;/Vo 


The ideal thrust loading coefficient of the propulsor is 
equal to: T 


Ст; = 1 е (5) 
5 уйла? 
zP 0 
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The thrust of the propulsor may be also obtained from 
the momentum conservation principle applied to the control 
surface S consisting of the side surface of the stream S, inlet 
cross-section of the stream S,, outlet cross-section of the stream 
S, and the surface containing the external forces acting on the 
fluid S, (Fig. 1): 


Pe . 48 = т Јрваз (& 


If the surface S, encompassing the external forces, is 
excluded from the surface S, then the axial force (thrust) T is 
equal to the rate of change of fluid momentum on the entire 
control surface and the force induced by the pressure field on 
the remaining part of S: 


T=- Граа = [РУСИ -mas + [pds © 
5р 5 S-Sp 
On the surface S, there is V- ii = — V, and on the surface S, 
there is V. ñ = V, + V,;. Consequently, the following relation 
between S, and S, may be established, making use of the 
continuity equation: 


MR} Vy = л (У, + V.) (8) 


In the case of an unbounded fluid domain the same pressure 
po acts on the entire surface S including surfaces S,, S, and 
S;, hence the second term in equation (7) is zero. In case of 
a bounded fluid domain (e.g. inside a tunnel) this condition 
does not hold. 

Taking into account that the side surface of the propulsor 
stream is in fact the stream surface, the scalar multiple У - ñ= 0. 
Finally, another formula for the axial force (or thrust) induced 
by an ideal propulsor is obtained: 


T = pVinR2m,4(l +) (9) 

Comparing formulae (9) and (4), the relation between 

the velocities induced at the propulsor and “far behind" the 
propulsor is obtained: 

2 dena 

M lud (10) 

It may be concluded from this relation that the induced 

velocity at the propulsor is equal half of the value of the 


induced velocity “far behind". Often this relation is accepted 
“а priori", which is not always correct, because it depends on 
the above assumptions, regarding an unbounded fluid domain 
and uniformity of velocity distribution over the propulsor 
surface. 

The thrust loading coefficient of the propulsor, related to 
the area of the propulsor stream cross-section “far behind” the 
propulsor is equal to: 


T 
Сту = 1 


— T (11) 
2 2 
5 РМ TR3 


=2m,,(1+m,,) 


The following relation stems from the principle of energy 
conservation (disregarding heat exchange): 
P,- [pv „505 + то ЈУ У -n)dS (12) 
S S 
The scalar multiple on the side surface of the propulsor 


stream V- їй = 0. The sum of integrals on the surfaces 5, and S, 
is also equal zero according to the continuity equation: 


[ром,45 x [роузаѕ = po( VIS, - V383)= 0 
51 53 
Hence the power (energy per second) lost in the stream is 
equal to: 
3 
1 2 V 
Py2-p [VV i)ds-pvjzR3| —9|- 2L = 
2 ч Vo 


1 
= pVonR; m3 (1 + т! + +m. 


The power loading coefficient of the propulsor, related to 
the area of the propulsor stream cross-section “far behind” the 
propulsor is equal to: 

Py; 1 

D = 
p3i ' = 2т „(1+ пој то (14) 
1 Мр 2 
2 P Vo TUN 3 

The same coefficient related to the propulsor surface area 

is equal to: 


Р {ү 
= Di = 
C»; T 1 ^ 2 = 2m + 2m.) (15) 
2 ОМ TR 2 
The efficiency of the propulsor is defined as the ratio of 
the power of thrust to the power lost at the cross-section “far 


behind” the propulsor (cf. Fig. 2): 
_ TiVo Crs 1 


-———— (16) 
Pp Съз 


| lem, 


1 
The efficiency of an ideal propulsor is often presented as 
a function of the thrust loading coefficients C,,; or Ст. 
Using the formula (11) a quadratic equation may be solved, 
leading to the following relation: 


which in turn may be used to develop the following formula for 
propulsor efficiency, depending on the thrust loading coefficient 
"far behind": 


(17) 


4 


(18) 


m3 10 


Fig. 2. The dependence of the ideal propulsor efficiency on the relative 
induced velocity "far behind" the propulsor 


In the case when the thrust loading coefficient at the 
propulsor surface is available, this relation takes a different 
form: 2 


Fig. 3 shows the dependence of the ideal propulsor 
efficiency on the thrust loading coefficients determined at the 
propulsor and “far behind". 

The above formulae refer to the optimum ideal propulsor, 
i.e. a propulsor operating in an unbounded fluid domain and 
having uniform distribution of thrust loading at its surface. 
When these assumptions are not fulfilled (this was the case of 
determination of the tunnel walls influence on the propulsor 
characteristics, as described in [3]), the induced velocity at 
the propulsor is no longer equal to half of the value of this 
velocity “far behind”. The deviation from this depends on the 
degree on non-uniformity of the thrust loading and on the ratio 
of the propulsor surface to the cross-section of the bounded 
fluid domain [3]. The vortex model of an ideal propulsor is 
composed of circular vortices distributed on the side surface 
of the propulsor stream behind the propulsor surface. 


(19) 


0 20 40 60 80 100 


2 120 Cu 140 


Fig. 3. The ideal propulsor efficiency as the function 
of thrust loading coefficients 


AN IDEAL FLUID BRAKE 


Physically an ideal fluid brake operates in the manner 
similar to an ideal propulsor, but now the pressure dipoles are 
directed in an opposite direction, i.e. against the direction of 
flow. Hence instead of an accelerating surface, there is now 
a decelerating surface and the axial force is generated in an 
opposite direction. The flow velocity is being reduced both 
in front and behind the brake surface (cf. Fig. 4). Similarly as 
before, the principles of momentum and energy conservation 
may be used for determination of the axial force and power. The 
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Fig. 4. The stream of an ideal fluid brake 


formulae are similar to those describing the propulsor, only the 
signs of induced velocity are changed and certain limitations 
in the values of these velocities are introduced. 

The axial force (now named the braking force instead 
of thrust) coefficient related to the stream cross-section “far 
behind" is now equal to: 


Сту = -2m,4(1- та) (20) 


and the coefficient of power loss in the stream, also related to 
the cross-section “far behind" is equal to: 


1 
Съз = 2m, «(1— пт) (21) 


Similarly as in the case of a propulsor, for the brake 
operating in an unbounded fluid domain and having uniform 
loading distribution over its surface, the velocity induced at the 
decelerating surface is equal to halfthe velocity induced at the 
cross-section “far behind" the brake. Hence the ratio between 
the areas of the brake surface and the stream cross-section “far 
behind" may be obtained from the continuity principle: 


ja 1 
nR? _ = Шз 
nRj 1-m,, 


Consequently, the coefficients of axial force and power 
loading, related to the brake area, have now the following 
form: 


1 
Cra o mim (22) 


2 
1 
Съ: = тај (23) 

The functions showing the dependence of the above 
coefficients on the induced velocity are plotted in Figs 5 and 6. 
The range of variation of the relative induced velocity “far 
behind" was limited to 1.0, because larger values would lead 
to physically unrealistic reversed flows. 

In case of a fluid brake it is inconvenient to discuss 
its "efficiency". The more appropriate parameter is its 
effectiveness. For example the ratio ofthe braking force power 
to the power loss in the brake stream may be proposed as the 
measure of the brake effectiveness: 
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=— = =i = (24) 


Fig. 5. Coefficients of kinetic energy loss in the brake stream depending 
on the relative induced velocity “far behind” m,; 
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Fig. 6. Coefficients of the braking force depending on the relative induced 
velocity “far behind” тз 


The value of the above parameter varies from 1.0 for my, = 
= 0.0 to 0.5 for my, = 1.0 . The above presented formulae 
contain all relations between kinematic and dynamic parameters 
of an ideal brake. 

It must be stressed here, that it is not correct to identify 
the value of C,,; determined by the formula (23) with the 
ideal efficiency coefficient for the axial wind or water 
turbines. This was suggested by Betz [1], who has used 
formula (23) for determination ofthe maximum ideal efficiency 


of а wind turbine (cf. Fig. 5). According to his result, quoted 
in almost all publications concerning wind turbines, this 
maximum value is equal to Сш, = 16/27 = 0.5926, which is 
obtained for: туу; = 2/3 = 0.6667. 

Contrary to the above, in experiments with real turbines 
having finite number of blades and operating in a real, viscous 
fluid (e.g air or water) significantly higher values of Cp, are 
obtained. This results from the fact that the operation of an 
ideal axial turbine should be determined not on the basis 
of the ideal fluid brake, but in the way similar to that 
presented below for an ideal screw propeller. The power of 
a turbine depends on the torque on its shaft and consequently, 
the circumferential induced velocity must be taken into 
account. In this case the principle of conservation of moment 
of momentum is more appropriate for determination of torque 
on the turbine shaft and the performance of the turbine should 
be made dependent on a parameter characterizing the moment 
of momentum. 

The energy conservation principle, taking into account all 
components of the induced velocities may also be used, but 
this way of solving the problem is much more difficult and less 
clear. Both ways lead to the same final result. 


AN IDEAL SCREW PROPELLER 


An ideal screw propeller model takes into account the 
rotational motion, therefore it must have a form of a circular 
disc. In order to determine the relations between the induced 
velocities and the forces on an ideal screw propeller, a more 
complicated vortex model should be used (cf. Fig. 7). 

This model is constructed of: 

- Straight line vortices distributed radially on the propeller 
disc, 

- helicoidal tip vortex lines distributed over the cylindrical 
side surface of the propeller stream behind the propeller, 

- Straight line vortex located along the propeller axis behind 

the propeller, having the total intensity I; 


The velocity field induced by such a vortex system at 
a certain distance behind the propeller (practically about one 
propeller diameter is enough) may be univocally determined by 
the parameters of the vortex system “Таг behind" the propeller. 
These parameters are: 
- radius of the cylinder on which the helicoidal vortices are 
distributed “far behind” R,, 


- total intensity of the vortices Г, 
- velocity of the undisturbed flow Уу, 
- pitch of the helicoidal vortex lines (ап, 


The volumetric mean value of the axial velocity induced 
inside the stream (r < R,) by the system of infinite helicoidal 
vortices has a constant value across the stream [4], equal to: 


Via =Г/(2 л R, tan.) (25) 


This value applies to a single helicoidal vortex line as well 
as to a system of Z (including Z = ©) heilcoidal lines having 
the same pitch equal to tanB,, as long as the total intensity of 
all vortex lines is equal to Г. At the same time the volumetric 
mean velocity induced outside the stream (r > R3) is equal 
zero [4]. 

The volumetric mean circumferential velocity induced 
inside the propeller stream (r < R,) by the infinite straight line 
vortex located at the propeller axis 1s equal to: 


Мев = Г/(2 лт) (26) 


The volumetric mean of the circumferential velocity 
induced by the system of tip vortices inside the stream (r « 
R,) is equal zero, while this value outside the stream (т > R;) 
18 determined by the formula (25) with inverted sign due to 
inverted direction of the tip vortices. Consequently, outside 
the stream the volumetric mean ofthe circumferential velocity 
induced by the entire vortex system of the propeller (axial 
vortex together with tip vortices) is equal zero. The volumetric 
mean of the radial component of the induced velocity “far 
behind" the propeller is equal zero both inside and outside the 
propeller stream. 

When the non-dimensional value of the intensity of vortices 
forming the propeller vortex system is introduced as: 


Ба = I/(4 x К, V) (27) 
the relative values (divided by the undisturbed flow velocity 


V) of the induced velocity components may be expressed by 
the following formulae: 


V/V, = My; = 2 by,/ tanp, 


МауМо = 2 буут = my, tanByr (29) 


When the velocity field “far behind" the propeller is 
defined in the above way, the relations between this field and 
propeller thrust T, torque Q and power loss may be developed 


Q8) 


——- 
* 


Vo*(1*m,;) 


Fig. 7. The stream of an ideal screw propeller with vortex lines 
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similarly as for an ideal propulsor, making use of the principles 
of conservation of momentum, moment of momentum and 
energy. 

It may seem that due to the constant value of the axial 
induced velocity across the stream, the value of propeller 
thrust may be determined analogically as in the case of an 
ideal propulsor. Unfortunately, the presence of the axial hub 
vortex significantly changes the pressure field inside the 
stream. Therefore, the thrust must be determined on the basis 
ofthe momentum conservation principle applied to the control 
surface S, which includes the side surface of the stream S,, 
inlet and outlet surfaces S, and S, together with the surface 
encompassing the external forces acting on the fluid Sp (cf. 
Fig. 7). 

[pv -ii)ds = - [pias (30) 
S S 

If the surface S, (representing the blades of the propeller) 
is excluded from the control surface S, then the thrust is equal 
to the rate of change of fluid momentum on the entire control 
surface and the force generated by the pressure field on the 
remaining part of S: 


T=- [pias = [pV(V-ü)dS [pads 1) 
Sp S S-Sp 

As there is V. її =- V, in cross-section | and there is V. ñ= 
= V, + Ма in cross-section 3, the relation between areas of 
cross-sections 1 and 3 may be determined using the principle of 
continuity offlow. However, it should be taken into account that 
on the axis of cross-section 3 there exists a vortex line, which 
should be treated as a vortex singularity (e.g. a Rankine vortex 

model), having cross-section area S,, and radius К. 


(S,-S,) (VotV,3) = S, V, (32) 


In the case of an ideal propulsor discussed above the same 
pressure p, acts on the entire control surface S (assuming 
unbounded fluid domain). In the case of ideal screw propeller 
the pressure in cross-section S, is different than p, and it is 
a function of radius, therefore an additional pressure term 
must appear in the corresponding equation. The side surface 
of the propeller stream is still a stream surface, where У ii = 0. 
Taking into account the above, the relation (31) now takes the 
following form: 


T= paR 3 МС + V,3)S;- Sw) + = ро)5 (33) 
53 


Substituting the relative induced velocity, this may be 
transformed into: 


T 2 pVenRAm,4(1- m4)(1-(R,, / R ))4 


34 
+ [E -poas di 
S3 


The first term of (34) is similar to the corresponding term 
in the formula (9) for an ideal propulsor. Calculation of the 
second term requires additional assumptions: 

- theentire vorticity in the flow is concentrated in the limited 
space around the tip and axial vortices; outside this limited 
space the flow is irrotational i.e. potential, 

- according to Joukovsky hypothesis the flow is steady in the 
system of co-ordinates rotating with the propeller. 


In sucha situation the following relation exists between the 
velocity field and the pressure field: 
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1/2V? - V e У, + p/p = const (35) 


where: — 
V.=@xR - the convective velocity (© is the angular 
propeller velocity), hence: 


p- po7 1/2 p( VÀ + VO + РСА Мр or. V) G6) 


Taking into account the equations (28) and (29) and 
reformulating the equation (36), the following equation for the 
thrust correction resulting from the pressure field in the screw 
propeller stream may be obtained: 


fe -pods = 
S3 


1 
=1/2 руту 2ле 5 Ја +tanB3/r? )r,dr,+ (37) 


TW 


1 
+ pV,m32nR; [a —г, tanB,/A, )г,аг, 
TW 
where: 
т, = r/R, 
= Vy / R; 
r, = R/R; is the normalized radius of the axial vortex 


The above integrals may be solved analytically and 
ultimately they lead to the formula for the thrust of an ideal 
screw propeller: 


Т =pVoaR3m,;((+m,,)(1—r )+ 


+ V, my; tan p; (2In(r,,) +1 –12 ) 


Now the thrust loading coefficient related to the propeller 
stream cross-section "far behind" is equal to: 


Стз = 1 Т т -2m,4(l* m X1- r )+ 
5 PVo mR; (39) 


(38) 


+ 4 түз tan В; 21n(r,) +1- n ) 

The power loading coefficient, related to the screw propeller 
stream cross-section “far behind”, may be determined from the 
principle of conservation of moment of momentum. In this case 
only the circumferential component of the induced velocity is 
important. This leads to the formula: 


Q,2 [p(POx У), (Ў -й)45+ |рРОха),45 (40) 


5-5 5-5 
where: 
PO - radius-vector from propeller axis 
ñ — unit normal vector at the control surface S-Sp 


An appropriate selection of the control surfaces enables 
elimination of some ofthe integrals. On the entire surface S-Sp 
the vector multiple (POxii), is zero, eliminating the pressure 
term. On the stream surfaces the scalar multiple V- ii is zero. 
On the cross-section S1 the following relations hold: 


V.n--V, 
(POxn), -[PÓx (-n-V)] 20 
(n), 2 -1 


and on the cross-section S, the following relations hold: 


V-A =V +V 
(РОхп), = PO- V, 
(m =1 


Taking the above into account, together with the continuity 
equation, the following 15 obtained: 


Qui =p „Мз (+ Vx3) dS 
833, 


After substituting the expressions (28) and (29) for the 
components of the induced velocity, the following equation is 
obtained for the power of an ideal screw propeller: 


Р, = 0, 00 = 
= pVonR 3 tanB3/A3 туз (1+шуз)(1— r2) 


where: 
) = У (о) 


(41) 


(42) 


The power loading coefficient, related to the area of the 
propeller stream cross-section “Таг behind", may be now 
determined as: 


1 
Cpi = P;/5PVonR3 = 


= 2tan B; /A3 my3 + туз)(1—гу) 
The efficiency of an ideal propeller is defined as the ratio 
of the power of thrust to the power loss in the propeller stream 
"far behind" the propeller: 


= ТУ = 
' P Di (44) 
=(1+1/2m,,, tan? Q In(r, )+1-12))A,/tanB, 


The formula (44) does not resemble the formula (16) for 
the efficiency on an ideal propulsor. However, it is similar 
to the formula used in the vortex theory of propellers for 
determination of the so called induced pitch for the propeller 
blades having optimum radial distribution of circulation: 


tgp = уп 


(43) 


(45) 
where: 
X= VyoR 


The formula (44) may take the form similar to formula (16) 
if the following definition of tanB, is used: 


вл = МЕМ 1+ 1/2; (46) 
У @R3-Vyy ИХ -1/2my; tanb; 
where: 


Vw and V „w — the components of the own velocity of the tip 
vortex. For an ideal propeller they are equal 
half of these components induced inside the 
propeller stream [4]. After some reformulations 
the following equation is obtained: 


1 
А, /їап = —— 
UE 1+1/2m,,, (1+ tanf2) 


After substitution of (47) into (45) the formula similar to 
that for the efficiency of an ideal propulsor is obtained: 


_ 1% _1+Y/2my, tanB5(2In(r,) +1= n) 
| Ppi 1+1/2my3(1+ (ап 87) 


(47) 


(48) 


The formula (48) differs from the formula (16) in two 
aspects: 
- the numerator is smaller than 1 
- the denominator includes an additional factor 


Both these differences cause reduction ofthe ideal propeller 
efficiency in comparison with an ideal propulsor. Therefore, the 
efficiency of an ideal screw propeller is always smaller than 
that of a corresponding ideal propulsor. Moreover, the ideal 
screw propeller efficiency depends additionally on the pitch of 
the helical vortex lines and on the diameter of the axial vortex 
singularity. Both formulae coincide when tanB, = 0. In this case 
the helical vortices of an ideal propeller are reduced to the ring 
vortices of an ideal propulsor. 

The radius of the axial vortex singularity requires further 
explanation. In case of an ideal screw propeller operating in 
water the cavitating vortex kernel may determine this radius. 
In case of an ideal screw propeller operating in gas (e.g. in air), 
there exists an exact relation between the pressure and volume 
(or density) based on a perfect gas model. In this case the axial 
vortex may be treated as a perfect whirlwind having a kernel of 
zero vorticity. The radius of such a kernel may be substituted 
into all above formulae in the place of rw. 


AN IDEAL AXIAL WIND TURBINE 


Similarly as in the case of an ideal propulsor and ideal 
brake, there exists a close analogy between an ideal screw 
propeller and an ideal axial wind turbine. Namely, they both 
may be described by the formulae characterizing the vortex 
system “far behind”: 

- radius of the cylindrical propeller stream “far behind” R,, 
- total intensity of the vortices Г, 

- undisturbed flow velocity V,, 

- pitch of the helical vortex lines tanf,. 


The vortex system of an ideal axial turbine, shown in Fig. 8, 
is composed of the same elements as the corresponding system 
for an ideal propeller: 

- Straight line vortices distributed radially on the propeller 
disc, 

- helicoidal tip vortex lines distributed over the cylindrical 
side surface of the turbine stream behind the turbine, 

- Straight line vortex located along the turbine axis behind 

the turbine, having the total intensity Г. 


Of course it must be remembered that the vorticity vectors 
for the propeller and for the turbine have opposite direction. 

In this case the velocity field is related to the parameters 
of the vortex system in the same way as in the ideal screw 
propeller. The relations (25) — (29) have the same form, but the 
parameters mv have opposite sign. The relations for the axial 
force T (38) and for the power Pi (42) have the same form, but 
the sign of m,, is changed: 


T = -pVznRm,,((l -m, (1 12) - 
(49) 


- lo my, tanB3(2In(r,,) +1-г„)) 


(50) 
=—pVonR; tan B;/A; ту my; )(1— гу) 


The definitions of all parameters аге the same as for an ideal 
screw propeller. Consequently, the formulae for the coefficients 
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Fig. 8. The stream of an ideal axial turbine with vortex lines 


of axial force loading Су; and power loading С, are similar to 
those for an ideal propeller: 


T 
Стз = TERE = -2т„5(( -m,)(1- 17) 
2 РУ TR; (51) 


-1/2 my, tan B2 (21n(r, ) -1— 12)) 


Сру =Р//2 руте з= 
(52) 
=—2 tanB;/A3 туз (1 - my3)1 - r) 
But now the formula for the pitch of the helical lines tanB, 
does change into: 


Vo + V. 1 -1/2m 
tan В» — 0 ху _ / V3 (53) 
from which the following relation may be developed: 
tanB,/A, = 1-1/2m,,(1+tan2p,) (54) 


Taking into account the formula (54) the relation (52) may 
be transformed into the form similar to (21): 


Са = -2 my (1-m,;)(1-1/2m,,(1+tan?,;))(1-1,7) (55) 


In this way it may be proved that in the system dependent 
on the axial component of the induced velocity m,, the limiting 
maximum values of the power coefficient С are even smaller 
than for the ideal fluid brake (assuming that ап, = 0 corresponds 
to the brake model and leads to the formula (21)). 

Taking into account that the axial velocity component 15 
only weakly related to the moment of momentum (i.e. to the 
turbine torque) and that the model experiments indicate values 
of C, higher than those obtained from the ideal fluid brake, the 
analysis ofthe turbine operation should be performed in another 
system of co-ordinates. For example, the coefficient b, more 
closely characterizes the moment of momentum in the flow. 

Considering the definition formula for m,, the formula for 
Cp; takes the form: 


Cj Ву2 ВУ ТЕЗ 
(56) 
= 463/03 @—2ЬууЛап 3;)(1— гу) 
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As in the case of a turbine the efficiency is not defined 
in the same way as in the case of propeller, the above 
formulae complete the relations between the dynamic and 
kinematic parameters. The radius of the vortex singularity 
may be determined in the above described way. Taking into 
consideration the continuity equation in the stream behind the 
turbine, the power loading coefficient at the turbine may be 
easily determined: 

Съд = 8,/0.5рмотЕ 5 = 
(57) 
--4by; /А2 (0 - 2by; /tanB;)l – та) 

This power loading coefficient is in fact the effectiveness of 
utilization of the energy of the oncoming flow by the turbine. 
The formula (56) enables further analysis of the influence of 
different parameters on the turbine efficiency described by 
Cy 
- the parameter (1 —r,,), characterizing the vortex singularity, 

has little influence on the turbine efficiency, 

- the parameter b,, = Г/4лУ,К, characterizes the circulation 
of the axial vortex and the moment of momentum in 
the stream. This is a very important parameter, which 
determines the torque and power of the turbine, 

- the parameter tanp, characterizes the pitch of the tip vortex 
lines at the turbine, which is difficult to be determined 
precisely, 

- the parameter А, characterizes the location of the turbine 
operating point on the co-ordinate of the advance 
coefficient J. 


It is difficult to determine the maximum attainable value of 
Cp; using the formula (57) and the above parameters. Clearly, it 
will be not equal to the Betz limiting value of Cpmax = 0.5926, 
determined for an ideal fluid brake. 

The structure of formula (57) implies that the maximum 
values of Cp are attained when: 

- the parameter 7, reaches minimum values, 
- the parameter Г describing circulation of the hub vortex, 
corresponding to b, reaches the highest values, 


Quantitative determination of these parameters is difficult. 
Certain idea of this aspect may be drawn from the analysis of 
the results of calculations for wind turbines with finite number 
of blades presented in the next section. 


Fig. 9. Experimental verification of the program СНАНҮР in application to wind turbines 


AN AXIAL WIND TURBINE WITH THE 
FINITE NUMBER OF BLADES 


Results of calculations showing that the effectiveness 
of a real axial wind turbines C, may exceed the upper limit 
defined by Betz have been obtained already in the 1990s, 
when the computer program for calculation of ship propeller 
characteristics CHAHYD was adapted for calculations of the 
axial wind turbines. 

The program CHAHYD, based on the vortex lifting surface 
theory, was very thoroughly verified experimentally in its 
application to ship propellers. Similar verification with respect 
to wind turbines was missing until 2010 [2]. Fig. 9 shows an 
example of such verification. In this figure the coefficients 
of power loading C, axial force K, torque K, obtained 
experimentally and computed for an 8-bladed turbine, are 
shown as functions of the advance coefficient J. In this figure 
the curve C, calculated for zero drag coefficient is also plotted. 
In this case the maximum value C, — 0.58 is already close to 
the limiting value determined by Betz. 

The experimentally verified results of calculations are used 
to formulate the following statement: 

The limiting value C max = 16/27 = 0.5926 as determined 
by Betz for an ideal brake may take significantly higher 
values for an ideal axial turbine. 


In the following figures the results of calculations of a 8- 
bladed wind turbine have been shown for four cases: 

- model wind turbine having pitch coefficient P/D = 1.0 
and viscous drag coefficient C, = 0.010 (these results are 
correlated with model experiments [2]), 

- the same model with P/D = 1.0 but with zero viscous 
drag, 

- themodel with pitch reduced to P/D — 0.8 and zero viscous 
drag (C; = 0.0), 

- themodel with pitch reduced to P/D — 0.6 and zero viscous 
drag (Cp = 0.0). 


As may be seen in Fig. 11 for the blade pitch P/D = 0.8 the 
С pmax Value already exceeds the Betz maximum of 0.5926. Still 
higher values of C max are obtained for a turbine with P/D = 0.6. 


If such values are obtained for a real wind turbine operating 


Fig. 10. Results of ране of the effectiveness с, 
for ап 8 bladed wind turbine 


Fig. 11. Calculated maximum values of C, obtained with viscous drag 
coefficients С, = 0.008 and С, = 0.0 as functions of the turbine pitch 
coefficient P/D 


in a viscous flow and having finite number of blades, then the 
corresponding values for an ideal wind turbine must be even 
higher. In Fig. 12 the results for wind turbines having different 
numbers of blades are presented. 

The quantitative assessment of terms in equation (57) 
Is difficult because the vortex model of a wind turbine with 
finite number of blades is markedly different from the ideal 
turbine with infinite number of blades. First of all, the fields 
of induced velocity are different and the limits of the turbine 
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С, = 0.008 
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Fig. 12. Calculated maximum values of C,, and C, (C, — 0.008) as functions 
of the turbine number of blades 


stream are not univocally determined. It is difficult to determine 

the parameters b, and my, for a turbine with finite number of 

blades. However, the above presented results of numerical 
calculations lead to the following recommendations for the 
wind turbine designers: 

- the design point of an axial wind turbine should be located 
at low advance coefficients, 

- the optimum pitch of a wind turbine blades lies near the 
value P/D = 0.6, 

- the computed effectiveness coefficient C, of a turbine 
increases with the number of blades, but for real turbines 
operating in viscous fluid and suffering from flow 
deformations near the hub, which grow with increasing 
number of blades, the number of blades should not be too 
high. 


CONCLUSION 


Analyzing the above results it may be concluded that the 
models of ideal propulsor and ideal screw propeller correlate 
well with the experimental data. The efficiency of real propellers 
is always lower than that ofthe ideal ones and the effectiveness 
coefficient (defined as the ratio of real propeller efficiency to 
the efficiency of an ideal one) characterizes correctly the losses 
connected with fluid viscosity and with the finite number of 
real propeller blades. 

The problem of turbines is significantly different, 
namely: 

- there is no correlation between the experimental results and 
the results of calculations, 
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- the analogy between a turbine and an ideal fluid brake, 
for which the axial induced velocity component m, 
is the basic parameter defining the performance, is 
a misunderstanding, 

- the turbine performance is characterized by its torque: 
the model of an ideal turbine reduces the values of C, as 
function of m, even further than the model of an ideal 
brake, but as it was mentioned earlier the torque 15 not well 
correlated with m, 

- when the results of model experiments are accepted as 
more trustworthy than calculations based on ideal models, 
the inadequacy of the Betz theorem for the determination 
ofthe limiting maximum value ofthe turbine effectiveness 
must be also accepted. 
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ABSTRACT 


Several different methods have been used to analyze the marine propellers. In this paper a boundary 

element method (BEM) for modelling steady hydrodynamic performance of a contra-rotating propeller 

(CRP) has been developed. The potential flow across the front and rear propellers have been considered. 

By calculating the induced flow, the interaction between two propellers is studied. This method has been 

applied to a typical CRP. The results are included in the form of the hydrodynamic performance coefficients, 

total thrust and torques. Numerical calculations indicate that the predicted hydrodynamic performance 
shows good conformity with experimental results. 


Keywords: boundary element method; contra-rotating propeller; high speed vehicle 


INTRODUCTION 


Because of increasing demands on higher load - carrying 
capacity and speed of the marine transport systems the marine 
research centres have concentrated their study on propeller 
devices that can meet this demand. However several problems 
arise when single propellers are used in large displacement ships 
or submarines to produce the required thrust. Amongst them 
are: increasing chance of propeller blade cavitation in heavy 
load conditions and vibrations due to unbalance in torque. The 
vibrations are transferred from the propeller to the ship body and 
cause instability of motion. The aforementioned problems must 
be considered by the designers specially in controlling rolling. 

A suitable method to overcome such problems is applying 
the contra-rotating propellers. Such propellers are composed 
oftwo propellers which are installed on two coaxial shafts and 
rotate in opposite directions. Despite the disadvantages such as 
complexity of design and heavy weight due to the propeller's 
motor, the contra-rotating propellers have many advantages 
over single propellers, including the following: 

a. generally, the rotational energy produced by the single 
propeller not only has no effect on the thrust but also 
causes loss of energy. The most important benefit of the 
contra-rotating propellers is their improved efficiency due 
to energy recovered by the back propeller. In fact, the back 
propeller recovers some amount of the rotational energy that 
is transferred to the water by the front propeller. Therefore 
using such propellers leads to a lower fuel consumption. 

b. the torque between the two propellers is balanced, that 15 
used to eliminate the reverse torque in devices in case when 


maintaining stability of motion is important. 

c. the load on the blades is lower and cavitation 15 delayed. 
The lower load also allows to reach a higher power output 
from CRP propellers of similar diameter compared with 
single ones. 

d. the produced noise is lower because ofthe uniformity ofthe 
downstream wake field of the contra-rotating propellers. 


Unfortunately there is no optimized design procedure for 
contra-rotating propellers and the full- scale test data are not 
published widely. 

There is no exact record of the history of contra-rotating 
propellers, but the first mention of them was given by 
Dollman and Perkins before 1800. Later, Ericsson performed 
experiments on model and real- scale propellers. Rota was 
the first who showed that the efficiency of the contra-rotating 
propellers is higher than the single propellers. He found that the 
reason for the increase in efficiency is the recovery of rotational 
energy transferred to water by the first propeller [1]. Several 
researchers performed numerical and experimental studies on 
the performance of the propellers and their parameters that 
affect it and their results were published in different papers 
some of which are listed by Cox et al. in [2]. Using lifting 
surface theory, Tsakonas et al. [3] computed the hydrodynamic 
forces applied on the contra-rotating propellers in stable and 
instable conditions. Yang [4, 5] used the lifting surface theory 
taking into account the interaction of the trailing wakes of the 
two propellers. Hoshino [6] also used the lifting surface theory 
to analyze the loads exerted onto the contra-rotating propeller 
shaft and has compared the results with the experimental 
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data gathered by him. In another research [7], Gu and Kinnas 
employed the Vortex-Lattice Method (VLM) and Finite Volume 
Method (FVM) to study the interaction of two component 
propellers in contra-rotating propellers and the propeller and 
duct (or nozzle) in ducted propellers, based on Euler solution 
of the entire flow field. Ghassemi [8] carried out numerical 
computations of the contra-rotating propellers for two types 
of ships, and recently [9] on the PBCF (Propeller Boss Cup 
Fin) effect to propeller performance. In his research, was 
also made a comparison between contra-rotating and single 
propellers for these ships. Koronowicz [10] described the use 
of computer software based on Vortex-Lattice Method. In the 
present paper his method 15 modified based on the boundary 
element method. 

As for experimental research the results of the experiments 
on a number of contra-rotating propellers were published 
by Hecker and McDonald [11]. Apart from evaluating the 
performance of each propeller, they also experimentally 
analyzed the effect of some parameters such as diameter of 
the front propeller, distance between propellers and pitch 
of the propeller. Van Manen and Oosterveld [12] performed 
experiments on a systematic series of contra-rotating propellers 
in which the front propeller had 4 blades and the rear propeller 
5 blades. The results of tests performed on a number of contra- 
rotating propellers, in which the loads applied on the contra- 
rotating propeller in instable non-uniform and uniform flow 
conditions were calculated, were also published by Miller 
[13]. Brizzolara et al. [14] presented a full numerical lifting 
line model with slipstream contraction effect to determine the 
chord pitch and camber distribution to respect the given margins 
of cavitation and strength. 

It is difficult to analyze a CRP in high rotational velocity 
field such as that of high-speed vehicle and there is no published 
data on such analysis using the boundary element method. In 
this paper, considering the physics of the problem we have 
obtained good results using the boundary element method. 
The results include the hydrodynamic performance of the rear 
and front propellers and their interaction. The computing of 
the induced velocity of the front propeller, regarded as the 
input velocity to the rear propeller, and the interaction of the 
two propellers are the most important problems which are 
considered in this paper. 


FUNDAMENTAL EQUATIONS 
Coordinate system 


We assume that the front and rear propellers rotate with the 
velocity o, and o, respectively, in non-viscous, incompressible 
and irrotational flow field. The uniform flow V, hits the 
propeller upstream in axial direction. 

We define the Cartesian coordinate system O-xyz (in which 
О 15 located at the centre of the propeller) such that x-axis 
is aligned with the propellers’ shaft and points downstream, 
z-axis is aligned with the generator line of the key blade, and 
y-axis is so located that a counter-clockwise coordinate system 
is formed (Fig. 1) 

In the cylindrical coordinate system, the angular position 
0 1s measured clockwise from z- axis. The coordinate r is 
measured from x- axis. Therefore, the Cartesian coordinate 
system is transferred to the cylindrical coordinate system by 
using the following equation is: 


х=х 
y = —rsinO (1) 
z-—-rcosO 
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Fig. 1. Coordinate system of СЕР 


in which: 


г= 4ју +z’ Q) 


Ө tan" (-y/z) 
Geometry of the propeller 


To construct the propeller, the geometry ofthe key blade of 
the rear propeller 1s formed in the Cartesian coordinate system 
and then transferred to the cylindrical coordinate system. Next, 
the Cartesian and cylindrical coordinates of other blades are 
obtained by using the coordinates ofthe key blade. For the front 
propeller the same procedure is repeated. As for the distance 
of the two propellers, the coordinates of the front propeller are 
displaced in axial direction. 

Usually, the cross-sections of the blade at different radii 
in 2D coordinate system are known. To define the geometry 
of the base cross-sections, the coordinate of each point P on 
the surface of the key blade is determined. For the clockwise 
(CW) and counter-clockwise (CCW) propeller the following 
equations could be used [15]: 


Xp=—[r tan(O, ) + 10, tan(B, )1+ (0.5 – x, )sin(B, ) 


+у cos(Bg) (3) 

| | 180[0.5 –х „)сов(Во) — Yup M! 

yp= rsinð t —— 
Tr 


2р= Г соја 


Tr 


In which, (+) applies to CCW propeller and (-) to CW one. 
The parameters В, 0, and 0, are: the geometric pitch angle, 
rake angle and skew angle, respectively. y,, is the distance 
from upper and lower surfaces to the chord line, obtained as 
follows: 


= у, + 
| y, ty, cosy "s 


YL =Y. Ур COSY 


where: 

(х„ у.) 
y(-0 * 5 te) 
Í 


— camber line coordinates, 
— the half thickness of the blade, 
— the slope of the camber line. 


For obtaining the coordinates of points on other blades, 
we use the following mapping in which ф is the phase angle 
of each blade: 


Xp 1 0 0 | xp 
У» |=|9 cosp sing | yp (5) 


Zp 0 sino cos@ | Zp 
Modelling the trailing vortex wake 


We assume that the surface of the wake, S,, is very small 
and there is not any discontinuity in pressure and flow on this 
surface. Discontinuity in potential values is allowable. The 
wake geometry should be modelled so as to satisfy the Kutta 
condition. As mentioned earlier, in order to simplify the solution 
procedure the thickness of the wake is assumed equal to zero. 
But the edge of each cross-section has some thickness and 
therefore the average of the coordinates of the top and bottom 
surface of the blade is used to determine the wake points on the 
edges of the blade. The total wake propagation in the direction 
of the blade axis is equal to the propagation of the propeller 
after two complete revolutions. Experimentally, this amount of 
propagation is sufficient to satisfy the Kutta condition. 

Changes in axial direction when the wake passes over the 
trailing edge of each section, are assumed negligible. It is also 
assumed that the propagation of each streamline of the wake is 
made on a constant-radius cylinder. In real situation, the radius 
of the trailing vortex wake (TVW) initially increases (transient 
region) and then decreases. 


The potential-based boundary element method 
(BEM) 
According to the made assumptions the disturbance 


potential satisfies the Laplace equation. Let’s consider the 
closed region V with the boundary S (Fig. 2). 


Fig. 2. Body mesh arrangement and its entering flow 


The boundary value problem for the potential equation 
could be expressed as Neumann problem as follows: 
V?$b-0 in V 
a (6) 


on S 


The kinematic boundary condition which should be satisfied 
on the boundary S,, is described by the following equation: 


A ўа where: V, 2 V, +@xF+Ving (7) 


Here, о 15 the propeller angular velocity of the propeller, V, 
is the advance velocity and Ving is the induced velocity of each 
propeller coming from the other propeller. 

To determine the vorticity the Kutta condition on the trailing 
edge of the blade should be satisfied. The Kutta condition is 
usually described as a finite velocity on the trailing edge. In 


other words: 
[Vore о (8) 


Considering the surface of the front and rear propellers 
and using Green theorem for the central point p on the surface 
of the front and rear propellers, S,, the potential problem is 
defined as: 


д 1 
27 (p) = [Ф dn, К(р,9) 


1 
— $, > dS + 9 
JM R(p.q) i 
0  l qs 
On, R(p.q) 


S, is the sum of the surface of the front propeller, Spp and 
that of the back propeller, Sga. Sy is composed of the surface 
ofthe wake of the front propeller, S,,,, and that of the wake of 
the back propeller, S,,,.The matrix form of Eq. (9) is [16]: 


[271 = Ом Фуа = 
= [Seas {n а UV hac {Аф ua 


Here, I is unit matrix, D, S and W are full matrices. ©, is 
determined by Eq. (7). The Kutta condition must be used for 
determining the unknowns Аф ofthe doublet strength for trailing 
sheet vortex surface. From Eq. (10) the velocity potential ф is 
calculated at each element. After the velocity potential over 
each element is determined, velocity and pressure distributions 
can be calculated directly by using Bernoulli's equation: 


v =Уф 


dS + 


+ JE лека) 


(10) 


Tr rr an 
P =0.5p ОМ м -¥,¥,) 
The pressure coefficient can be expressed as: 
р 
C, mL (12) 
0.5 pV, 


Thrust and torque of the propellers with two components 
of pressure and friction are expressed as follows: 


N 
T=Z)> Pn AS, -T, 


4 Fric 
Е (13) 
Q= Ly Pr -n,,-¥, JAS, + не 


i=l 


where Ту and О; are the frictional component of thrust and 
torque of the propeller, respectively. They are determined by 
using empirical formulae proposed by ITTC. The total thrust 
and total torque produced by the CRP is given by: 


Те 
о= 0 |-[0,] 


(14) 
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Hydrodynamic coefficients 


The hydrodynamic coefficients of each propeller can be 
obtained. For the front propeller we have: 


T Qr 
Кү = = ; K of bo: NN 
pn, D; pn Dj б 
пе = Krf J J= Va 
| Кој 27 n,D, 


In which n,1s the rotational velocity, Dis the diameter, T, is 
the thrust and 0,15 the torque produced by the front propeller. 
For the rear propeller we have: 


T Q 
Ка = а > Ко = Е 
Ta pn, ^ D,* Qa pn, D,? 
Kg J (16) 
C Ko 27 


In which n, 15 the rotational velocity, D, is the diameter, T, is 
the thrust and Q, is the torque produced by the rear propeller. 

Finally, as the thrust and torque of the rear and front 
propellers and performance coefficients of each of them are 
known the performance coefficients ofthe CRP can be obtained 
as follows: 


Ky = Ку +Къ 


Ko = Коу +К о 2 
da 
Ko 2n 
Induced velocities 


From Green's theorem, as applied to the potential field in 
Eq. (9), we can alternatively construct the velocity field by 
taking the gradient of the perturbation velocity: 


Fig. 3. Inflow velocities of CRP 


ATV ing = 118 $(q)V n, POR 
1 
= M (18) 
БЕ 
+ Д АФУ г. em 
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The induced velocities due to the front propeller, directed 
to the aft one, is defined by Ум, and vice-versa - by V, qr; 
Therefore, the inflow velocity to the aft propeller is obtained 


as follows (Fig. 3): 
Vi = V4 + Op XT + Vinda) + Уша(АЕ) (19) 
V2 = Va + XT + Vina AF)  VindFA) 
The action of the rear propeller onto the front one is very 
low compared with the action of the front propeller onto rear 
one, hence it can be neglected in some cases. 


NUMERICAL RESULTS 
Hydrodynamic performance 


In this section, the results of analysis of a typical CRP are 
presented. The analyzed propeller is of a 3686F-36874 type. 
The front propeller (3686F) is a CCW propeller and the rear 
propeller (3687A) is a CW one. The geometrical parameters 
of the propellers are given in [13]. Main parameters of both 
the propellers are presented in Tab. 1. An image of the whole 
СЕР 15 shown in Fig. 4. The trailing vortex wake (TVW) is 
shown in Fig. 5. 


Tab. 1. Main dimensions of the CRP 


Parameters 3686F 3687A 
Blade number | 4 O 
Diameter [m] 0.3052 0.2991 


P/D at 0.7R 
EAR 
Rotation 


0.303 
Anticlockwise 
NACA66 a = 0.8 


Clockwise 
NACA66 a= 0.8 


Blade section 


Run ofthe hydrodynamic coefficients of the front propeller 
(3686F) is shown in Fig. 6. and for the rear propeller - in 
Fig. 7. The coefficients for both the propellers (1.e. whole 
CRP), compared with experimental results, are shown in 
Fig. 8. The numerical results show good agreement with the 
experimental data. 

In 2003, Gu and Kinnas [6] analyzed this type of propeller 
(3686F-3687A) using coupling vortex lattice and finite 
volume methods. In Fig. 9 1s presented a comparison made 
between results of the present study, Gu and Kinnas research 
and experimental one. It can be observed that the results of 
the present study are within a small error range, hence rather 
acceptable. 


Fig. 4. 3686F -36874 СЕР model 
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Fig. 6. Hydrodynamic coefficients of the front propeller 
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Fig. 7. Hydrodynamic coefficients ој the rear propeller 
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Fig. 8. Hydrodynamic coefficients of the 3686F-3687A CRP model 
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Fig. 9. Comparison of the numerical results (of the present study and 
Kinnas [6]) and experimental data for the considered CRP 


Calculations of the thrust and torque 


Values of the ratio of the thrust produced by the front 
propeller and the rear one related to the total thrust at the relative 
radius r/R = 0.7 are obtained, respectively, as follows: 


TR, = 81е _ 47 279 
di Thrust total 
Thrust 
ТЕ да = A = 52 7396 
Thrust total 


It can be seen that the thrust produced by the rear propeller 
is higher. This may be due to the effect of the front — propeller 
- induced wake on the rear propeller. 

Fig. 10 and 11 illustrate the thrust-RPM, torque-RPM, 
total thrust-RPM and total torque-RPM relations, respectively. 
Vectors of the torque produced by the front and rear propellers 
are directed against each other and thus the total torque are 
much lower than the individual torque of each propeller. 
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Fig. 10. Individual thrust and torque of each propeller versus propeller 
rotational speed (RPM) 
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Fig. 11. Total thrust and torque of the whole CRP versus propeller 
rotational speed (RPM) 


Pressure distributions on blades and induced 
velocities 


The coefficient of contour pressure distribution on the 
propeller blade surfaces at J — 1.19 are shown in Fig. 12 
and 13. Red colour shows the high pressure face side and blue 
one indicates the low pressure back side. The induced velocities 
are very important because the interactions between two 
propellers are caused by them. The axial induced wake velocity 
between two propellers and that behind the rear propeller is 
shown in Fig. 14 and 15, respectively. 


Fig. 13. Pressure distributions on the back (upper or left?) and face (lower 
or right?) surface of the rear propeller 
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Fig. 15. Axial induced wake velocity behind the rear propeller (X/R = 0.34) 


Fig. 12. Pressure distributions on the back (upper or left?) and face (lower 
or right?) surface of the front propeller 
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CONCLUSIONS 


With regard to the numerical results of the 3686F-3687A 

CRP calculations the following conclusions can be offered: 

• As for the hydrodynamic characteristics of the front and 
rear propellers, it could be mentioned that the thrust and 
torque hydrodynamic coefficients for the rear propeller are 
higher than those for the front one. This is due to the effect 
of the front propeller wake. 

* The hydrodynamic performance parameters of the CRP were 
obtained by superposition of the parameters of individual 
propellers. This is consistent with the experimental 
results. 

* The CRP minimizes the torque produced by the propeller; 
hence it is recommended to apply it to devices for which 
elimination ofroll movement is important (e.g. torpedoes). 
Also eliminating the torque causes less vibration to be 
transmitted to the ship body and thus increases the stability 
of motion of the ship. 

* As for the pressure distribution on the rear and front 
propellers, it could be concluded that in a CRP the load is 
divided between front and rear propellers. Each propeller 
has its own contribution in producing the total thrust. 
Therefore using this system results in smaller loads on 
each propeller and consequently more uniform pressure 
distribution on the blades. Therefore in contra- rotating 
propellers chance of cavitation appears lower. 


NOMENCLATURE 

D — propeller diameter 

D; — fore propeller diameter 

D, — aft propeller diameter 

EAR - expanded area ratio 

Z — number of blades 

J — advance velocity ratio 

P/D — propeller pitch 

Кү — thrust coefficient of the fore propeller 
К — thrust coefficient of the aft propeller 
Ку — thrust coefficient of the propeller 

Ко — torque coefficient of the propeller 
Koa — torque coefficient of the propeller 

Ko — torque coefficient of the propeller 

P — pressure 

C, — pressure coefficient 

Sw — trailing wake surface 

Sp — blade surface 

T — total thrust 

T; — fore propeller thrust 

T, — aft propeller thrust 

Tas — frictional component of thrust 

Q — total torque 

Q; — fore propeller torque 

Q, — aft propeller torque 

Qui — frictional component of toque 

Va — advance velocity 

Vs V, — downstream velocities of the propellers 
V, — inflow velocity to the propeller 

у, — tangential induced velocity 

Vind — induced velocity 

Vaga) — induced velocity by fore propeller to aft propeller 
Vnaar) — induced velocity by aft propeller to fore propeller 
D, S, W — potential matrix coefficients 

n, — fore propeller rotating 

n, — aft propeller rotating 


n — normal vector of the surface 


Yu — distance from upper and lower surfaces to the chord line 
(X, у.) — camber line coordinates 

te» — half thickness of the blade 

Or — angular velocity of fore propeller 
Q, — angular velocity of aft propeller 
Вс — geometric pitch angle 

0 — rotating angle of blade 

Өк — rake angle 

05 — skew angle 

Ф — phase angle of each blade 

у — camber line slope 


ф — velocity potential 

Da — normal derivative of potential 

Аф — difference potential at trailing edge 
in — inflow velocity potential 

AS — surface area of each element 

n — propeller efficiency 

un — fore propeller efficiency 

Na — aft propeller efficiency 
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ABSTRACT 


Trim control mechanisms such as interceptors and trim flaps have been widely used in recent years in high- 
speed crafts for ride and trim control. In spite of their extensive application, a few studies investigating 
the impact of interceptors on planing craft performance, have been published. In the present study, the 
impact of interceptors on planing crafts hydrodynamic quality is investigated through application of an 
experimental method. Two scaled-down models of high-speed planing mono-hull and catamaran are tested 
with and without interceptors in calm water at different heights of the interceptors to investigate the effect 
of interceptors on drag reduction of the models. The first one is a scaled-down model of 11 m planing 
mono-hull boat and the test was conducted at the towing tank of Sharif University of Technology, Iran. 
The second one is a scaled-down model of 18 m planing catamaran boat and the test was conducted at 
the towing tank of Krylov Shipbuilding Research Institute (KSRI), Russia. The experimental results show 
a remarkable drag reduction of up to 15% for mono-hull model and up to 12% for catamaran model over 
the wide speed range of the models. 


Keywords: interceptors; drag reduction; hydrodynamic quality; planing catamaran; mono-hull 


INTRODUCTION 


Interceptors have been used for trim and ride control 
of planing and semi-planing boats since the 1970s and the 
application of these devices has steadily increased over the 
years. Interceptors are nowadays also used for steering and 
motion control of high speed crafts. An interceptor is a thin plate 
or blade, protruded from the hull normal to the flow direction 
causing a stagnating flow region near the blade. The resultant 
pressure acts on the hull bottom creating the effects such as the 
trim moment, which adds lift and finally controls the attitude 
of the craft. The interceptors are normally installed at the boat 
transom and/or amidships or at a distance ahead of transom 
depending on the boat loading, position of the centre of gravity 
and initial trim. Hydraulic or electric actuators adjust the height 
of the interceptors and the blades can be retracted completely 
when they are not needed. The most common application of 
interceptors is to serve as a trim control device that can help 
to reduce the overall resistance of planing vessel. There is also 
another way to decrease the resistance of a planing vessel by 
means of interceptors, i.e. through installing the interceptor 
amidships. This results in the reduction of the wetted area, 
both by increasing the dynamic lift and creating an air cavity 
behind the interceptors. The magnitude ofthe lift created by the 
interceptors normally depends on their height and flow velocity. 
However, there are other parameters possible to have an 


influence on the effectiveness of the interceptors. Fig. 1 shows 
the effect ofthe arrangement ofthe interceptors on modification 
of lift and wetted area of the boat by combination of added lift 
and reduced wetted surface brought about by the interceptor's 
effect and whereby the resistance of the boat can be reduced 
and attitude of the vessel be controlled respectively. 
However, besides the extensive application of the 
interceptors, there have been a relatively few published 
studies either on investigating the impact of interceptors on 
hull resistance or on a systematic study on the performance of 
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Fig. 1. Illustration of physical phenomena associated with interceptors 
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these devices. Thus, the present study aims at improving the 
understanding of the impact of interceptors on hull resistance 
by reviewing the published literature sources, examining the 
impact of the interceptors on hydrodynamic quality of planing 
mono-hull and catamaran through an experimental study, 
and to develop further insight into the procedures by which 
interceptors can cause reduction of resistance. 


LITERATURE REVIEW 


The effect of stern wedge on ship powering performance of 
planing and semi-planing vessels has been studied by several 
researchers [1]. Bannikov, Bannikova et al. experimentally 
investigated the transom interceptor influence on the 
hydrodynamic characteristics of planing surfaces with 0.3 m 
beam, dead rise angle of 0, 15 and 30 deg at 4-9 m/s speed, at 
different heights of the interceptor, and derived an empirical 
formula for the added lift and drag of planing surfaces [2, 3]. 
By using several planing surfaces and flat plates'experimental 
results, Savitsky devised a computational procedure to evaluate 
the performance of a planing hull. The foregoing procedure is 
used in calculating the running trim and resistance at a given 
running speed as well as the load and location of the centre 
of gravity [4]. 

Like interceptors, the trim flap has also been widely used on 
semi-planing and planing vessels and there is a large amount 
of published data on its performance. A thorough investigation 
was conducted by Brown as "an experimental and theoretical 
study of planing surfaces with trim flaps" [5]. Furthermore, the 
effect of the trim flaps was implemented to the planing boat 
model by Savitsky and Brown [6]. 

Tsai et al. carried out a model test on a 1/20 scale model 
of 20 m patrol boat at designed speed of 40 knots and a 1/10 
scale model of a 29.5 m patrol boat at designed speed of 32 
knots, including interceptors and with and without stern flaps. 
The application of even small interceptors (h/L of the order of 
0.196) changed the trim and as a result reduced the resistance 
considerably (up to 17% at interceptor height of 0.5, 1, 1.5 mm 
in model scale) [7, 8]. 

The Krylov Shipbuilding Research Institute also performed 
several tank tests on catamaran and mono-hull models to 
explore the effectiveness of interceptors on resistance reduction 
and motion control of high speed crafts [9, 10]. 

A new formulation has been developed by Dawson 
and Blount as the basis to predict the equivalent lift of 
interceptors [11]. 

Brizzolara carried out one of the most detailed studies 
on interceptor's hydrodynamics. He used a CFD method to 
study the local flow around a 2D interceptor fitted to a flat 
boundary, which represented the hydrodynamic flow around 
the bottom of a ship and was the most idealized model of 
interceptors ever used. This model allowed for viscous effects 
and through limited series of CFD-based simulations provided 
the relationships between the hydrodynamic forces and the 
interceptor's main parameters [13]. Successively, Molini and 
Brizzolara introduced a simple potential flow model for the 
predicting of pressure and lift forces in front of the interceptors 
[14]. D.Villa and Brizzolara put forward a new study, which 
changes context and approach considerably in comparison with 
the previous ones. In fact, CFD-based hydrodynamic models 
consider the complete description of planing hull form with 
stern appendages. Of course, not only interceptors but also 
stern flaps have been considered with the scope to compare 
relative performance; the CFD solution with RANES solver 
were used to simulate the flow around planing hulls and stern 
appendages more precisely [15]. 
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In a series of M.Sc.-theses elaborated in the Department of 
Marine Technology at Norwegian University of Science and 
Technology (NTNU), Norway, their authors tried to establish 
a database for performance data on lift and drag of interceptors 
installed at the transom of prismatic planing hulls. Sverre 
Steen summarized the experimental results suggesting various 
approaches to empirical formulas for the coefficients of added 
lift and drag due to interceptors [16]. 

Alterskjar performed the tank test of the model of 26 m 
catamaran with stern and mid-interceptors in calm water and 
waves to investigate whether the application of stern and 
mid-interceptors could help reducing the calm water resistance 
and if the good results observed in calm water are also present 
in seaway [17]. Based on his finding, Allema has reported 
a successful attempt in reducing the required power for a cruise 
ship [18]. 

Fridman has presented the most thorough research 
material available on interceptors used as a means to ventilate 
ship's hull bottom. He also reported the research results 
on a system consisting of both under-bottom-mounted and 
stern-mounted interceptors on cutter-class vessels. Fridman 
also suggested theoretical models to predict lift and drag due 
to the interceptors [19]. 

Alexander H. Day et al. carried out a comprehensive 
experimental tank test program on a scaled model of a yacht 
hull to examine the impact of interceptors. The results indicate 
a significant reduction in calm water resistance over a wide 
speed range with gains of 10 — 18 % in the speed range between 
8 and 20 knots, accompanied by a reduced trim and sinkage. 
The performance of the interceptors was also compared with 
Gurney flap, and in addition the effect of interceptors on trim 
change was compared with the longitudinally moving ballast, 
and as observed, the benefits appear to be largely sustained in 
small waves [21]. 

The literature review of the available results on interceptor 
and its application indicates that in spite of the popularity 
of interceptors, there have been a relatively low number of 
published studies comprehensively exploring the effect of 
interceptors on hull resistance, and apparently, the foregoing 
surveys were not specifically aimed at high speed planing 
mono-hull and catamaran resistance. Therefore, the aim of the 
present study is to provide a method to evaluate the effect of 
interceptors on high speed planing mono-hull and catamaran 
and fill the knowledge gap that currently exists, by carrying 
out the tank test of the scaled-down planing mono-hull and 
catamaran models with interceptors at full planing speed. 


SELECTING INTERCEPTOR SIZE AND 
DIMENSIONS 


Interceptor height 


Several researchers have studied the optimal height 
of interceptors to maximize the interceptor effectiveness 
in terms of improving hydrodynamic quality of planing 
vessels. Most of them agree that the selection of interceptor 
height depends mostly upon the boundary layer thickness 
at interceptor location. The idea widely adopted is that the 
interceptors should be contained entirely within the boundary 
layer, so that the interceptor height is partly a fraction of the 
boundary layer thickness. Thus, as the initial approximation, 
the boundary layer near the transom (for aft interceptor 
location) can be considered in the same way as that over the 
flat plate having the same length as the hull. The thickness 
for the boundary layer in the turbulent flow can be calculated 
from Eq. (1), [22]: 


8(х)= 0.37 xRe, 4 (1) 


Several manufacturers producing interceptor systems have 
put forward various suggestions for the interceptor height. 
Humphree suggested the interceptor height of up to 50 mm 
for power boats between 18 and 45 m, [25], while other 
manufacturers recommended a height extended to 75 mm for 
heavier vessels between 18 and 60 m in length. 

Brizzolara used an interceptor at the maximum height of 200 
mm on the steering interceptor of STENA HSS-1500 vessel of 
127 m overall length and 40 kn speed [13]. Alexander H.Day 
et al. surveyed the literature and plotted the found interceptor 
height values against the Reynolds number. The mentioned author 
compared then the results with the boundary layer thickness values 
obtained from Eq.(1), [22]. The relevant figure is re-plotted in Fig. 
2 with the additional data, showing non-dimensional interceptor 
height in relation to the boat length. Since tank tests usually cover 
a wide range of speeds, only several values are plotted for each 
study, which represent different interceptor heights and different 
speeds in order to put the scatter in context. Because of the fact 
that in the case of planing vessels the wetted length can vary 
substantially with speed, the static wetted length is used. 
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Fig. 2. Experimental data on non-dimensional height of interceptors 


As depicted in Fig. 2, the selected interceptor's heights 
become much smaller proportionally to the boundary layer 
thickness. Based on the above-mentioned suggestion and 
as shown in Fig. 2, in the present experimental study the 
interceptor was tested at model-scale height of 2 to 3 mm for 
mono-hull model and 1to 4mm for catamaran model. Therefore, 
the non-dimensional interceptor height varies from 0.4 96 to 
0.6 % of the model waterline length for mono-hull model, 
and from 0.03 % to 0.13 % of the model waterline length for 
catamaran model. As Fig. 2 shows, at corresponding Reynolds 
numbers, the selected interceptor heights are located well within 
the range of the boundary layer thickness (from 13.6 % to 
23.7 % of the boundary layer thickness for mono-hull model 
and from 2.8 % to 13.5 % of the boundary layer thickness for 
catamaran model, respectively). 


Interceptor spans and positions 


In addition to selection of interceptor height, the other 
important parameters for the case of the location of both stern 
and mid-interceptors are: span and longitudinal position of 
interceptors. As for the span, the full beam of the mono-hull 
and the half beam of each twin-hull (130 mm) of the catamaran 
model are used as interceptor's span. For both the models 
transom is considered as the location of the stern interceptor 
and the longitudinal position of the mid-interceptor for the 
catamaran model is the distance of about 3 = 5 % of waterline 
length ahead of the model’s centre of gravity, [12]. 


CATAMARAN MODEL EXPERIMENTAL 
PROCEDURE 


Model and model tests 


The 18 m planing catamaran selected for the tank-test 
investigation, has V-shaped cross sections and hard chine 
twin-hulls. In the current study, all tank tests were performed 
in the KSRI high-speed towing tank ofthe dimensions: 650 m 
(L) x 14 m (W) x 6 m (D). Its carriage is able to travel at the 
speed ofup to 16m/s (from 0.01 m/s up to 16m/s with 0.01 m/s 
step) with the speed measurement accuracy of +/-0.001 m/s. 
During the tests, the model was fixed against sway and yaw, 
while free to heave, trim and roll. The 1/6 scale model was 
built of wood and a special lightweight foam and formed 
with the use of CNC milling machine followed by special 
paint spraying procedures to give perfect finishing to the hull 
surface. The model had special windows in the bottom of each 
of twin-hulls for monitoring wetted surface boundary behind 
the mid-interceptor. The experimental procedure in general 
conformed to the ITTC standard procedures for model making 
and resistance testing [23, 24]. During the towing tests the 
following parameters were measured: speed of the model, total 
resistance of the model, dynamic trim and dynamic sinkage at 
the point of towing strut attachment. In addition, the boundary 
of hull-wetted area was also recorded at each test speed. The 
model’s principal dimensions are given in Tab.1. 


Тађ.1. Models principal dimensions 


Model Boat Main particulars 


3m 18m Length overall 
0.70 m 4.2m Beam (total) 
65 Kg 144t Displacement 

11.55 m/s 55 knots 


(0.32 + 0.38) x L 
measured from 
transom 


20 deg 
1.56m 
1.08 m 


Design speed 


(0.32 + 0.38) x L 
measured from 
transom 


20 deg 
0.26 m 
0.180 m 


Гер, (centre-of- 
gravity location ) 


Dead-rise at Lcg 


Twin-hull beam 


Tunnel beam 


Test matrix 


As Fig. 1 shows, in the case of considering location of both 
stern and mid-interceptors, the relative position of interceptors 
and the boat centre of gravity will define the state of dynamic 
stability ofthe vessel, and before implementation of any control 
mechanism the vessel itself should be dynamically stable in 
terms of porpoise motion etc. Then in the present study, the 
optimum position of the model’s centre of gravity was selected 
on the basis ofthe bare hull performance, and next the tank test 
was performed for the selected position ofthe centre of gravity 
and various interceptors’ height values. 


Performance of bare hull at different 
centre-of-gravity positions 


The model of 65 kg in weight was tested in bare hull 
condition for three relative positions of centre of gravity 
corresponding respectively to 0.32L, 0.35L and 0.38L measured 
from transom ahead, over the speed range corresponding to full- 
scale speeds up to 60 knots. Fig. 3 and 4 show the variation of 
relative resistance and trim angle along with the displacement 
Froude number (FNA), respectively. Based on the results 
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a number of observations can be made: the variation of the Leg 
of the model has remarkable effects on the relative resistance 
and trim angle over the wide range of Froude number values 


corresponding to the planing speeds. 
0.3 - T" — T — = = 


0.25 


0.2 


0.15 


~@Leg=.32xL -dHLlcg-.35xL --Lcg-.38xL 


Resistance/Weight Ratio 


0 06 12 18 24 3 36 42 48 54 6 6.6 
Displacement Froude Number ЕМА 
Fig. 3. Variation of relative resistance with Froude number 
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Fig. 4. Variation of running trim angle with Froude number 


Tab. 2 shows the variation of the relative resistance and 
trim angle along with FNA. As Tab. 2 shows, at the Froude 
number corresponding to the first hump (ЕМА = 1.5 ~ 1.8), the 
shift of the Lcg up to 6 %L to the bow will give the resistance 
reduction up to 11 = 14%. 

As Fig. 3 shows the second hump is observed at FNA = 4. 
Atthis FNA value which corresponds to 38 kn speed in the full 
scale (lower than cruising speed), the maximum resistance is 
observed at the Leg = 0.38L measured from transom. From 
FNA = 4 up to ЕМА = 5.82 (corresponding to the full scale 
design speed), the resistance will be decreased by shifting the 
Leg to the bow. From the dynamic stability point of view, as 
observed during the tank test, at the relative position of the 
centre of gravity, corresponding to 0.32L, the model suffered 
porpoising instability at ЕМА > 6 + 6.2. Based on the above 
given results, the relative longitudinal position of the centre 
of gravity corresponding to 0.38L measured from transom 
was then selected for execution of tank test of models with 
interceptors. 


Tab. 2. Variation of resistance and trim with ЕМА 


Relative resistance 
reduction ratio 


0.32L | 0.35L | 0.38L 


0.88 
(11.594) 
0.86 
(1496) 
0.95 
(5%) 


Relative dynamic 
trim reduction 


0.32L | 0.35L | 0.381. 


ЕМА = 1.0 


First Hump 
ЕМА = 1.5 «1.8 


ЕМА = 3.0 


Test-setup for interceptor and actuator system 


Resistance test in calm water is carried out not only with the 
fixed-height interceptors but also the selected relative position 
of centre of gravity (Lcg = 0.38L) at the model displacement 
of 65 kg. 
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The stern interceptors of the height h,,,= 1 mm, ћ = 2mm 
and h,,,=4 mm and the mid-interceptors of the height h;,,= 1 mm 
and ћа = 2 mm, were tested. The stern interceptors were 
mounted on the transom board. The distance between the stern 
interceptors and the mid-interceptors was 1.3 m. 

The interceptor pull-rod was inserted through an opening 
in the model hull and connected to a lever, with hinge axle. 
Another lever mounted to the same axle was connected to 
a lanyard attached to a foundation on the model's deck. By 
turning the lanyard the lever axle is caused to rotate and finally 
the interceptors is brought about to slide vertically along the 
hull side. a pointer mounted to the axle indicates the extended 
height of the interceptors. 

Prior to each test run, the interceptor's height settings were 
calibrated within 0 and 6 mm. This calibration established the 
relation between the position of the axle-mounted pointer and 
the interceptor's extended heights measured vertically. 

The interceptors were made of 4 mm thick Al-alloy 
plates. Lower edges of the plates followed the relevant 
portion of the hull lines. In order to reduce the friction of 
the movement of the interceptors, the bottom step surfaces 
in the way of interceptor's attachments were flush-covered 
with the vinyl plates. Each interceptor was secured to the 
respective step by two bolts. The interceptor's shapes and 
the attachment details are shown in Fig. 5 and Fig. 6 depicts 
the model hull lines. 


Fig. 6. Catamaran model hull lines drawing 


Results of calm-water tests with fixed-height 
interceptors 


During the calm-water tests of the model having the 
displacement A — 65 kg and centre-of-gravity location 
Leg = 0.38L, the interceptors were located individually either 
at the middle or transom section, with various interceptor 
height settings simultaneously at the both sections. Testing 
a series of the fixed-height settings was sufficient for finding 
the optimum interceptor height. The test was performed for the 
complete speed range corresponding to the range of FNA up 
to 6.0. Measurements made during the tests were as follows: 
model speed, hydrodynamic resistance, and trim at towing 
strut attachment point. Besides, observers visually marked the 
foremost boundary of the portside twin-hulls wetted keel and 
at the chine. The information were necessary to predict the 
wetted area for subsequent full-scale extrapolations by using 
the Froude method. 

Plots in Fig. 7 through 12 demonstrate the interceptor 
height and their configuration effects in relation to the relative 
resistance (R/A), the trim and the rise of the centre of gravity, 
at FNA = 1.8 and FNA = 5.6, corresponding to the hump and 
design speed of the full-scale vessel, respectively. 
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Fig. 7. Variation of the relative resistance with the interceptor height at 
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Fig. 8. Variation of the rise of the centre of gravity with the interceptor 
height at hump speed 
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Fig. 9. Variation of the running trim with the interceptor height at hump 
speed 
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Fig. 10. Variation of the relative resistance with the interceptor height at the 
design speed 
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Fig. 11. Variation of the dynamic trim angle with the interceptor height at 
the design speed 
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Fig. 12. Variation of the rise of the centre of gravity with the interceptor 
height at the design speed 


A number of observations can be made based on the 
obtained results: 

When the model runs with the hump speed (FNA = 1.8), 
extending only the mid-interceptors by 2 mm (2, 0), 
h/L=0.06 %, causes a sizeable resistance growth and the model 
was unable to ride faster than with FNA = 3.6 because of the 
onset of intensive oscillations. 

Simultaneous 2 mm extension, (h/L = 0.06 %), of both the 
mid- and stern- interceptors (2, 2) brings a certain decrease in 
the resistance only at FNA = 1.8, whereas at higher speeds the 
resistance abruptly starts to grow. 

Extending only the stern-interceptors by 2 mm (0, 2), 
h/L= 0.06 %, reduces the resistance significantly compared to 
riding without interceptors at low speeds. Further increasing the 
interceptor height to 4 mm (0, 4), h/L = 0.13 %, had virtually no 
effect on the resistance at the hump speed. In accordance with 
this interceptor height setting (0, 4), by increasing the speed 
the model trim angle and rise of centre of gravity decreased 
drastically. 

During higher-speed tests (ЕМА = 4.5 + 6.0), interceptor 
height settings were varied within the range of 1 + 2 mm. 

Extending the interceptor height to 1 mm (h/L = 0.03 %) 
for the mid-interceptor only (1, 0), similar to the extending of 
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interceptors at the both sections to the same height (1, 1), results 
in increasing resistance at higher speeds. 

Extending the interceptor height to 1mm (h/L = 0.03 %) at 
the transom only (0, 1), reduces the resistance at the moderate 
Froude number values, FNA = 4.5 + 5.6. 

Based on the analysis of the obtained results, the most 
efficient measure in terms of reducing the resistance at the 
hump and higher speeds was to extend the stern interceptors to 
2 mm (h/L — 0.06 %). In some cases, the resistance reductions 
were as much as by 10 — 12 %, accompanied by noticeable 
reduction oftrim angle. For the higher speed values (FNA > 6), 
extending the interceptors resulted in a sizeable resistance 
growth and the model became unstable (started performing 
intensive oscillating motions), therefore it is not advisable 
to use the fixed interceptors when the model rides with such 
high speeds. 


MONO-HULL MODEL EXPERIMENTAL 
PROCEDURE 


The 11 m planing mono-hull designed by the authors and 
selected for the tank-test investigations, has V-shaped cross 
sections and hard chine. In the present study, all tank tests 
were performed in the marine laboratory towing tank, Sharif 
University of Technology, with the dimensions of: 25 m (L) x 
2.5 m (W) x 1.25 m (D). Its carriage is able to travel at speeds 
up to 6 m/s (from 0.01 m/s up to 6 m/s with 0.5 m/s step) with 
the speed measurement accuracy of +/- 0.001 m/s. During 
the tests, the model was fixed against sway and yaw, and free 
to heave, trim and roll. The 1/15 scale model was made of 
wood and a special lightweight foam and formed by means 
of CNC milling machine, followed by special paint spraying 
procedures to ensure perfect finishing to the hull surfaces. The 
experimental procedure in general conformed to the ITTC 
standards, procedures for model making and resistance testing 
[23, 24]. During the towing tests, the following parameters were 
measured: speed of the model, total resistance of the model, 
dynamic trim and dynamic sinkage at the point of towing strut 
attachment. In addition, the boundary of hull-wetted area was 
also recorded at each test speed. Principal dimensions of the 
model are given in Tab. 3. 


Tab. 3. Model principal dimensions 


Main particulars 


0.73m 


11.0m 


Length overall 


0.16m 


2.40 m 


Beam 


2.02 kg 


6.831 


Displacement 


5.31 m/s 


40 knots 


Max speed 


0.35 x Lk measured 
from transom 


15 deg at transom, 
18 deg at centre of 


gravity 


0.35 х Lk mea- 
sured from tran- 
som 


15 deg at transom, 
18 deg at centre of 


gravity 


Longitudinal posi- 
tion of centre of 
gravity 


Dead-rise angle 


Model initial trim — 1.25 deg 
Model keel wetted length( Lk) = 0.6 m 
Model chine wetted length (Lc) = 0.4 m 

Model draught at transom — 0.05 m 


Fig. 13 shows the model hull lines and Fig. 14 presents the 


model and related instruments in the towing tank. 
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Fig. 14. Model and related instruments in towing tank 


Test matrix 


The model having 2 kg weight and one relative position 
of centre of gravity corresponding to 0.35L measured from 
transom, was tested in bare hull condition over a speed range 
corresponding to full-scale speeds of up to 45 knots. Fig. 
15 and 16 show the variation of relative resistance and trim 
angle along with displacement Froude number, respectively. 
The stern interceptors of the height ћ, = 2 mm and h,,, = 
— 3 mm were tested. The stern interceptors were mounted 
to the transom board. The interceptors were made of 4 mm 
thick special plastic plates attached to transom by means of 
a special paste. 

Details of shape of the interceptor are shown in Fig. 17. 
Resistance test in calm water was carried out not only with 
fixed-height interceptors but also with the selected relative 
position of centre of gravity (Leg = 0.35L), and the model 
displacement of 2 kg. 
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Fig. 15. Variation of the relative resistance with Froude number 
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Fig. 16. Variation of dynamic trim angle with Froude number 


Fig. 17. Details of shape of the interceptor 


Results of calm water tests with fixed-height 
interceptors 


During the calm-water tests ofthe model ofthe displacement 
А = 2 kg and the centre-of-gravity position Leg = 0.35L, the 
interceptors were located at the transom section, with various 
height settings. The test was performed for the complete 
speed range corresponding to the FNA values of up to 5.3. 
Measurements made during the tests were as follows: model 
speed, hydrodynamic resistance and trim at towing strut 
attachment point. Besides, observers visually marked the 
foremost boundary of the portside hull wetted keel and at the 
chine. The information was necessary to estimate the wetted 
area for subsequent full-scale extrapolations by using the 
Froude method. Plots in Fig. 18 through 21 demonstrate the 
interceptor height and its configuration effects in relation to the 
relative resistance (R/A) and trim at different Froude number 
values. 
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Fig. 18. Resistance test results for two interceptor s height 
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Fig. 19. Drag reduction effect due to 2 and 3mm interceptor 
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Fig. 20. Dynamic trim test results for one interceptor height 
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Fig. 21. Variation of the relative resistance and trim with Froude number 


A number of observations can be made based on the 
obtained results: 

Based on Fig. 15 it is found that the first resistance hump 
occurs at FNA = 1.56 which corresponds to the maximum 
running trim angle of 6 deg (Fig. 16), as Fig. 16 shows the 
running trim at the design speed of 40 kn (FNA = 4.8) is equal 
to 3.7 deg. 

Fig. 18 depicts the relative resistance of the model with two 
different interceptor height values. Fig. 19 illustrates the drag 
reduction effect compared to the bare hull test results. Fig. 20 
shows the effect of the interceptor on running trim. As Fig. 18 
shows, the effect of the interceptor on relative resistance of 
the model is different in terms of Froude number. As shown 
in Fig. 18 three different areas correspond to various Froude 
number ranges; each of them can be considered in evaluating 
the interceptor effectiveness. When the model runs at hump 
speed (area A, ЕМА = 0 + 1.5), extending the interceptors by 
2 mm, h/L = 0.4 %, and 3 mm, h/L = 0.6 % causes a small 
resistance growth resulted from low lift/drag ratio of the 
interceptor at such low speed values. When the model runs 
under transition to planing speed (area B, 1.5 < FNA < 3), 
extending the interceptors significantly reduces the resistance 
compared to riding without the interceptors. Extending the 
interceptor by 2 mm, h/L = 0.4 %, reduces the resistance up 
to 5.2 %, and extending it by 3 mm, h/L = 0.6 %, lowers the 
resistance by 10.5 % at FNA = 1.56. When the model runs 


POLISH MARITIME RESEARCH, No 2/2013 27 


at fully planning speed (area C, 3 « FNA), extending the 
interceptors still reduces the resistance of the model, but the 
sign of the resistance reduction in relation to the interceptor 
height, is changed at some Froude number depending on the 
interceptor height. Except that Froude number, extending the 
interceptor causes a sizeable resistance growth. Tab. 4 shows 
the percentage of resistance reduction for different interceptor 
height values. 


Tab. 4. Resistance reduction at different interceptor height extension values 


Resistance reduction 
percentage 
at h/L=0.4% 


3.3+5.5% 
7 * 11% 
11% 


Resistance reduc- 
tion percentage 
at h/L = 0.6% 


8+ 1% 
15 + 17% 
17 + 19% 


Range of FnA 


Extending the interceptor at higher speeds causes a sizeable 
reduction in running trim, and results in a noticeable resistance 
growth. As Fig. 19 shows, the increase in resistance occurs at 
FNA = 4.7 for interceptor extension of 3 mm, (h/L = 6 %) and 
ЕМА = 5.2 for interceptor extension of 2 mm, (h/L = 4 94), 
respectively. Then at higher speeds the rate of resistance 
growth 18 increased by a rise in the interceptor height. At 
the design speed of 40 knots, the resistance reduction was as 
much as 12 % for 2 mm interceptor extension, whereas 3 mm 
interceptor extension results in increasing resistance. Fig. 20 
shows that extending the interceptor by 2 mm, (h/1— 4 %) leads 
to a remarkable reduction in running trim. In fact, the trim 
reduction is about 10 % when the model rides at the hump speed 
and about 30 % when the model rides at the design speed. 


CONCLUSIONS 


The following conclusions can be obtained from the model- 
test results: 

1. Based on the present analysis and achieved results, it 15 
proved that the interceptors can be applied as a highly 
effective device to reduce the resistance of high- speed 
planing crafts, and in addition their performance can be 
improved by appropriate varying, in an integrated manner, 
Leg and size of the interceptors. 

2. Use of the stern- and mid- interceptors can give the 
remarkable reduction in the resistance and running trim. 

3. It is clear from the present study that for higher speeds 
extending the interceptors involves a sizable resistance 
growth. 

4. As for the catamaran model, the extension of the stern 
interceptor by 2 mm results in some cases in the resistance 
reduction of up to 10 + 12 96. 

5. As for the mono-hull model, the extension of the stern 
interceptor by 2 mm results in the resistance reduction of 
up to 7+ 11 % in some cases, and for its extension by 3 mm 
the resultant reduction in resistance was up to 15 + 19 %. 

6. Performing the experimental studies and developing 
the mathematical models for evaluating effectiveness of 
automatically controlled interceptor in calm water and in 
waves could be valuable future research tasks aimed at 
making the application of interceptors more effective in 
whole range of speed and loading conditions. 


NOMENCLATURE 


ӧ(х) — Turbulent boundary layer thickness [mm] 
R — Reynolds number at relative length position 


ex 
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h;,/L, h/L — non-dimensional interceptor height/length ratio 


int 


d/L — non-dimensional turbulent boundary layer thickness 


R/W — resistance/weight ratio 
FNA — displacement Froude number (== 
ү8У С 
V — ship speed [m/s] 
V — volume of displacement [m?] 
g — acceleration of gravity in [m/s?] 
А — boat or model displacement 
Leg — longitudinal position of centre of gravity as a fraction 
of boat length, measured from transom. 
ac h — interceptor height. 
L — ship or model length. 
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Wake-blade interaction in steam turbine stages 


Jerzy $wirydczuk, Prof. 
Institute of Fluid-Flow Machinery, Gdansk, Poland 


ABSTRACT 


The article discusses the phenomenon of stator Wake/Rotor cascade (W/R) interaction 
in a steam turbine stage, and the ability to capture it in turbine stage design calculations 
making use of standard numerical codes. Firstly, the W/R interaction is analysed by 
comparing its real, experimentally recorded course with the numerical results obtained 
using vortex theory models and methods. This part of the analysis ends with formulating 
a conclusion about stochastic nature of the W/R interaction and indicating its reason, which 
is the vortex structure of the stator wake. Next, a question is discussed whether and how 


this stochastic nature of the examined phenomenon can be taken into account in calculations of Reynolds 

Averaged Navier-Stokes (RANS) equations. Differences are indicated between the uniform pattern of 

the stator wake obtained using a RANS code and the vortex structure of the real wake. It is concluded, 

however, that despite these differences the RANS results correctly reflect the time-averaged course of the 

real W/R interaction, and the process of averaging the flow parameters on the sliding plane between stator 

and rotor calculation areas can be treated as sort of "numerical averaging" of different realisations of 
the W/R interaction. 


Key words: turbine passage; secondary vortices; stochastic behaviour; CFD calculations 


SECONDARY VORTEX STRUCTURES IN 
A TURBINE ROW 


Steam turbines are widely used in shipbuilding industry 
as the main engines on ships which need huge driving power 
(large super-tankers, for instance), or the nature of their 
operation is characteristic for long-distance sailing without 
calling at a harbour (ice-breakers or marine submarines with 
nuclear propulsion systems). Permanent competition between 
turbine producers forces them to offer high-quality and high- 
efficiency products, the working characteristics of which are 
to be presented to the buyer at an earliest possible stage of 
negotiations. These characteristics are to be calculated very 
precisely, as delivering the product which does not meet 
parameters named in the contract would result in fines paid by 
the producer. On the other hand, when the sold turbine reveals 
the real power output or efficiency much higher that that 
declared in the offer, it means an extra bonus for the buyer, or, 
looking from the other side, loss for the seller who could have 
sold it for a higher price. That is why steam turbine designers 
attempt to calculate basic turbine parameters very accurately, 
taking into account as many factors as possible. One of these 
factors which were given special attention in recent years 15 
the effect of interaction of the secondary vortices forming 
in turbine blade stator and rotor passages on the level of the 
generated stage loss. 
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The flow through a steam turbine stage is extremely complex, 
Fig. 1 [1]. Even in high pressure stages of relatively simple 
radial geometry in which the flowing superheated steam has 
physical properties close to those of the perfect gas, a number 
of secondary vortices, both forming in stator or rotor cascades 
and approaching from upstream rows, can be observed. Looking 
from the cascade front, two pairs of horseshoe vortices 1 and 2 
are formed at the leading edge of each blade, one pair near the 
hub endwall and one near the tip endwall. Then, inside the rotor 
passage these vortices mix, partially or entirely, with the passage 
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Fig. 1. Secondary flow structure in turbine stage: 1, 2) pressure and suction 
side horseshoe vortices, 3) passage vortex, 4, 5) suction and pressure side 
corner vortices, 6) trailing shed vortex 


vortices 3 generated by the passage cross flows. In the corners 
between the blade hubs, or tips, and the passage endwalls, the 
action of the passage vortices leads to the formation of corner 
vortices 4 and 5. The rear parts of the stator and rotor blades 
are the sources of wakes which shed from their trailing edges 
and move downstream into next turbine parts. The interaction 
of the passage vortex with the blade trailing edge 15 believed 
to be a source of a so-called trailing shed vortex 6, having the 
rotation opposite to that represented by the passage vortex [2]. 
Flow separations, occasionally observed at rotor blade surfaces, 
can lead to the creation of large-scale vortices of various 
orientations. And last but not least, leakage flows over the blade 
tips in the rotor cascades usually generate characteristic vortex 
structures the intensity of which depends on the adopted type 
of sealing systems. 

The scale and direction in which the interaction of all above 
vortices affects turbine performance parameters is not clear. 
Generally it is assumed that neglecting this phenomenon leads 
to underestimation of stage losses [ 3], but opinions can also 
be found that in some cases the presence of vortex structures 
improves the stage performance by periodical sweeping out 
the stagnation areas. 

Although extremely complex, the abovementioned vortex 
interaction in a steam turbine stage can be successfully analysed 
using Computational Fluid Dynamics (CFD) codes which 
solve the Reynolds Averaged Navier-Stokes (RANS) equations 
in computational domains of arbitrary geometry. It should 
be stressed, however, that these CFD calculations require 
extremely fine grids to take into account the widest possible 
range of small-scale vortices. A commonly accepted approach 
here is to perform the calculations in the geometry consisting 
of one stator and one rotor passage in combination with special 
procedures taking into account the relative rotor/stator motion 
when passing the flow field parameters between these two 
computational domains. The assumption which is believed to 
justify this approach 15 that all unsteady phenomena take the 
identical course in each stator and each rotor passage. 

However, there is one exception to this assumption, which 
Is the interaction of the wakes shedding from the stator blade 
trailing edges with the downstream rotor blade cascades, the W/R 
interaction. A detailed analysis of this phenomenon апа its effect 
on the turbine stage design process is the basic subject of this 
article, which synthesises the experience gained by the author 
in this area over the recent decade. Due to the lack of reliable 
experimental data on W/R interaction coming from marine 
steam turbines, in order to keep the geometrical and physical 
conditions of the analysis within realistic regimes the geometry 
of the cascade used in the analysis was assumed identical with 
the midspan cross-section of a real HP impulse turbine rotor in 
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Fig. 2. Negative jet concept of wake model 


operation in one of inland power plants in Poland. The selected 
parameters of the flow through the rotor cascade were also 
taken from the measurements performed in this turbine. The 
conclusions resulting from the presented analysis are valid for 
both inland and marine applications of steam turbines. 


WAKE MODELLING 


In the past, the wakes shedding from turbine blades were 
modelled in various ways depending on a current state of 
knowledge on the phenomenon. In one of the first models 
proposed by Meyer in 1958 [4] the wake was represented by 
a regular band of velocity deficit, a so-called "negative jet", 
immersed in the otherwise uniform flow, Fig. 2. The stator 
wakes were transported with the main flow and chopped into 
linear segments by the rotor blades. Within the rotor passage 
the wake segments, keeping their linear shapes, generated 
secondary velocity fields responsible for fluid transport from 
the pressure side to the suction side of the passage. 

Due to a very simplified geometry assumed for both the 
stator wake and the rotor blades, the Meyer's model could not 
take into account any possible deformation of the stator wake 
on its way through the rotor. However, as was noticed by Smith 
[5], the velocity difference between the rotor blade pressure and 
suction sides in a real turbine row makes the opposite ends of 
the wake segment move with different velocities, which leads 
to stretching and shearing ofthe entire segment. Consequently, 
when the wake leaves the rotor, its structure resembles a saw 
rather than the initial linear and continuous shape. 

Taking into account the wake deformation was only 
possible after the development of computer technology and the 
appearance of CFD codes which solve the fluid flow equations 
in two and three dimensions. Initial attempts to capture the wake 
deformation based on the potential flow model with a continuous 
vortex layer of changing intensity shed from the blade trailing 
edge, which was an adaptation ofthe classical Prandtl's concept 
of unsteady flow over an airfoil [6] to the cascade geometry 
[7, 8]. Once shed to the flow, successive wake segments moved 
with local flow velocities through the rotor cascade, thus making 
a further contribution to unsteady flow effects in this area, Fig. 3. 
Instantaneous intensities of the vortices, being the response to 
changing flow conditions in the vicinity of blade trailing edges, 
were calculated from a set of integral equations. 
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Fig. 3. Vortex wake model 


It is noteworthy that despite a discrete representation, 
the wakes in these models were still interpreted classically 
as continuous bands of fluid in otherwise inviscid and 
incompressible flow, with the velocity deficit changing 
uniformly along the wake axis. 

More advanced studies of the stator wake behaviour 
in the rotor cascade were possible using specialised 3D 
codes developed to solve RANS equations for turbulent, 
compressible, and viscous flows through fluid-flow machines. 
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In the calculations performed with the aid of those codes the 
stator wakes were obtained as part of the solution and their 
characteristics did not have to be assumed a priori. The RANS 
codes enabled studying numerous aspects of W/R interaction, 
including wake deformation [9], profile losses [10], unsteady 
heat transfer [11], boundary layer transition [12], pressure fields 
and blade load fluctuations [13], etc. Experimental verification 
of the obtained results was, generally, good — at the level 
determined by the accuracy and resolution of the measuring 
devices used in those times. 

The initial calculations performed with the aid of3D RANS 
codes made use of grids of relatively moderate resolution. 
As aresult, the generated wakes had the forms of regular 
layers of increased entropy and decreased flow velocity, 
without visible traces of inner structure, Fig. 4. However, the 
experiments oriented on taking more detailed insight into the 
stator wakes revealed clear frequencies in the signals recorded 
in the wake area, which testified to the presence of regular 
fine vortex structures in there [14, 15]. More recent numerical 
computations performed for isolated stators using sufficiently 
fine grids and time steps also produced wake patterns in which 
two rows of well developed vortices having the structure 
similar to that of the von Karman vortex street was observed 
[16 = 19]. In this context, a question should be raised about the 
consequences brought by the presence of active vortices in the 
wake to conventional CFD analyses making use, as a rule, of 
RANS codes as basic tools. 
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Fig. 4. Stator wake deformed inside the rotor passage — NS calculations 


VORTEX STRUCTURE OF THE WAKE AND 
ITS CONSEQUENCES 


Analysing numerically the development and interaction 
of the vortex stator wake with the rotor cascade is much more 
difficult that in case of all other secondary structures described 
in section 1. Particular wake segments chopped off by the rotor 
blades do not have, in general, the same characteristics, i.e. do 
not consist of vortices of identical strengths and positions at the 
time of chopping. And any differences in the initial structure of 
the chopped-off wake segment result in inevitable differences in 
its further motion through the rotor passage. Unlike the passive 
band-like wakes observed in conventional CFD calculations, 
the sample of which is shown in Fig. 4, active vortices interact 
both with each other and with the boundaries, which makes 
their trajectories differ from those estimated on the basis of 
pure main flow convection. 

The initial arrangement ofthe vortices in the wake segment 
atthe time of chopping-off by the rotor blade may be described 
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using a so-called phase shift, ø, which defines the range of 
displacement of the von Karman vortex structure along the 
wake axis. The definition of the phase shift is shown in Fig. 
5 where it 18 interpreted as the ratio of the distance s between 
corresponding vortices in two different von Karman wake 
vortex segments a and b to the distance / between two adjacent 
vortices of the same rotation in one segment. 
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Fig. 5. Phase shift in the stator wake vortex arrangement 


The above aspect of the W/R interaction cannot be 
effectively analysed using a NS code. Due to possible 
differences in phase shifts of the particular chopped-off wake 
segments we cannot assume any longer that their behaviour 
inside each rotor passage is the same. Consequently, to 
capture the complete W/R interaction and its effects we should 
generate a computational grid covering all stator and rotor 
blade passages. And what is even more challenging, this grid 
should be fine enough to detect and study the behaviour of such 
small structures as the wake vortices. Keeping in mind that 
the number of passages in the stator or rotor row can exceed 
a hundred in a real turbine, the creation of a sufficiently fine 
grid in those circumstances is beyond capabilities of even most 
powerful computing centres. 

That is why studying vortex related aspects of the W/R 
interaction needs a flow model which would take into account 
the differences between vortex arrangements in individual 
wake segments. The model applied in the present study makes 
use of the Vortex Dynamics Theory (VDT) which defines 
principles of motion of vortices in multiply connected regions 
[20]. Along with the ability to model arbitrary phase shifts 
of the wake vortex configurations, VDT provides several 
invariants of vortex motion which can be used as additional 
tools to increase the accuracy of flow calculations. The here 
applied VDT model bases on the 2D potential flow through an 
infinite cascade of blades, with the wake represented by vortex 
singularities. Unlike the earlier works making use ofthe vortex 
model [7, 8], in the present study each singularity represents 
an individual coherent vortex. 


FLOW MODEL 


A procedure which is frequently applied in studying 
unsteady vortex phenomena bases on generating, step by step, 
a series of instantaneous vortex configurations to illustrate the 
course of the phenomenon. In this procedure each previous 
configuration is used as the input data for calculating the next 
one. Therefore of crucial significance for the accuracy of the 
final result is the accuracy with which transient configurations 
are calculated. 

In most typical vortex model applications in blade flow 
analyses the shape of the blade is modelled using a sequence 


of discrete vortices representing the bound vorticity distributed 
along the blade contour - see, for instance, Basu and Hancock 
for a single profile [21], or Lienhart [7], and Krammer [8] 
for a cascade. The above blade description can be a source 
of relatively large inaccuracy in evaluating instantaneous 
velocities of free vortices, especially those moving close to 
the blade contour. To reduce this inaccuracy, the present study 
makes use of conformal mapping [22] and the circle theorem 
based inversion technique [23] which move all singularities 
used for defining the blade shape off the boundaries, thus 
creating a continuous and regular velocity field in the entire 
flow area. A detailed procedure of constructing the stream 
function for a system of blades in various configurations can 
be found in Refs. [24 + 28]. The accuracy of the solutions 
obtained using the conformal mapping is believed to converge 
asymptotically to the accuracy of the analytical solutions. 
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Fig. 6. Infinite blade cascade mapped into a circular cascade of circles 


The mapping function used in the present analysis maps 
a K-element segment of an infinite linear cascade of blades into 
a circular cascade of K circles, Fig. 6. A general formula for this 
function has the form of an infinite fraction series [22]: 
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The number K of the blades composing the transformed 
cascade segment can be chosen arbitrarily. The resultant 
length of the cascade segment is T = tK, where t is the cascade 
pitch. In numerical applications the infinite sum of mapping 
coefficients is reduced to a limited number, N, the choice 
of which compromises between the required accuracy and 
the available computation time. Unknown parameters in 
Eq. (1) are the centre of the first circle б, the circle radius a, 
equal for all circles, and a series of complex coefficients 
C,. For the analysis performed in the rotor reference frame 
these parameters are to be calculated once in the entire flow 
calculation process. 

In the canonical plane, G the parallel flow transformed 
from the rotor cascade plane, z, is represented by the source 
Q, and the vortex Г, located at the origin S. The strengths of 
these two singularities are related with the assumed inlet flow 
velocity, w,, and inlet flow angle, p, in the physical plane and 
are calculated as: 
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for the vortex and: 
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for the source. 

A general formula for the complex flow velocity in the 
canonical plane б is also given in the form of a fraction 
series [22]: 
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Here Оу „аге unknown complex coefficients calculated from 
the system of equations expressing the condition of zero normal 
velocity along the entire contour of each blade, and zero tangent 
velocity at circle points corresponding to the blade trailing 
edges. These coefficients represent the reaction ofthe circles to 
both steady and unsteady flow components. S is the number of 
pairs of vortices in a single stator wake, and Г, is the strength 
ofthe positive vortex in the stator wake, a significant parameter 
in the present study. E and С, | are instantaneous locations of 
the positive and negative vortex, respectively, in the I-th stator 
wake. W is the number of vortices in one rotor wake and G, 15 
an instantaneous location of the rotor wake vortex generated in 
the w-th step downstream of the k-th blade. АП those parameters 
are to be calculated at each step of computation. 

The flow velocities calculated in the canonical plane are 
transformed to the physical plane using the formula [22]: 
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where 47/4 is calculated from Eq. (1). 


The velocity v; of the i-th vortex situated at point G in the 
canonical plane (here 1 stands for all indices of both stator 
and rotor vortices) can be calculated from Eq.(4) using the 
Helmholtz and Kirchhoff's principle that the vortex does not 
act upon itself [23]: 


_ = у il; 
vi(5;) = v(G;) - lim — — (6) 
сэс б - G6; 
The vortex velocities are transformed from the canonical 
plane to the physical plane using the formula [23] : 
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The derivatives in Eq. (7) are calculated from Eq. (1). 

In the VD Theory, the strengths and initial positions of the 
wake vortices are to be assumed a priori based on the available 
information on the examined flow, including in particular the 
scale of the velocity deficit generated by the wake and the 
width of the wake. Unfortunately, the literature on the subject 
gives little and rather incomplete data on the inner structure 
and strength of the vortex wakes in any turbine stages. In 
this situation, in order to obtain the wake characteristics 
which could be considered reliable and representative for the 
examined turbine flow, in the reported study the initial wake 
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parameters were reconstructed from preliminary 3D RANS 
calculations done using a specialised code FlowER for the 
stage being the prototype for the rotor geometry examined in 
the present study. 


STEADY FLOW TESTS 


Before the basic part of the analysis, the applicability of 
the VDT flow model to the examined phenomenon of vortex 
interaction inside the turbine rotor passage was checked on 
a series of comparison tests performed using the results of the 
FIowER calculation as the reference. 

The first test consisted in comparing the velocity distributions 
inside the rotor passage obtained using the potential flow model, 
without vortices, and FlowER. The input data for the potential 
flow calculations were the rotor inlet velocity w;, and angle 
В, computed by FlowER and averaged at the rotor inlet mid- 
span. The results of the test are shown in Fig. 7 as velocity 
modulus and angle distributions in selected passage sections 
located at: x/c = - 14.85%, 25.73%, 66.52% and 101.68% from 
the rotor blade leading edge. Bars in the diagrams represent 
the FlowER data, while the points located close to their tops 
stand for the results obtained using the potential flow model. 
The presented data well match each other, except for the rotor 
blade boundary layer near the blade suction side and the rotor 
wake where the viscous effect are dominating. 


X4/c = -14.85% хз/с = 66.52% 
Xolc = 25.73% X4/c = 101.6896 
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Fig. 7. Distributions of steady flow velocity modulus (lefi) and angle (right) 
in rotor passage 
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Fig. 8. Stator wake deformation in the rotor passage: 3D CFD midspan 
solution (lefi) and potential flow model (right) 


The second test was oriented on comparing instantaneous 
images of wakes computed by FlowER with those obtained 
using the potential model. A sample FlowER wake is shown in 
Fig. 8 (left). The corresponding potential model of the passive 
wake, Г, = 0, was generated by placing two parallel lines of 
fluid particles upstream of the rotor and then allowing them 
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to move with the flow through the rotor. Two instantaneous 
locations of these pairs of lines, selected so as to correspond to 
the wake locations generated by FlowER, are shown in Fig. 8 
(right). The agreement between the both curves is surprisingly 
good, bearing in mind complex structure ofthe flow in a turbine 
stage. These two tests testify to good credibility ofthe potential 
model to determine accurately the velocity fields inside the 
rotor passage, and thus to create a realistic pattern ofthe wake 
motion through it. 


THE W/R INTERACTION 


Following the standard procedure mentioned in the section 
entitled Flow model the course of the stator wake/rotor cascade 
interaction was studied by calculating a sequence of instantaneous 
flow patterns showing current distributions of stator wake 
vortices inside the rotor passage. The input parameters for these 
calculations were: the inlet flow velocity w,, the inlet flow angle 
В, and the strengths of the stator vortices Г,. The main variable 
parameter ofthe study was the phase shift factor, «, defined in Fig. 
5 (the basic distribution of the wake vortices, corresponding to 
Фф = 0.0, is labelled “a” in this figure). 

Figure 9 shows a selection of these flow patterns covering 
one period of wake passing. The flow is from left to right. Stator 
wake vortices are represented in the diagrams by small circles, 
green and red for the positive — counter-clockwise and negative 
— clockwise rotation, respectively. Black circles downstream 
of the rotor blades mark initial segments of the rotor wake. 
The diagrams also include two sequences of small black dots 
representing hypothetical positions of the stator wake vortices 
in case they moved passively with the potential flow velocity, 
like the fluid particles in Fig. 8, right. 

In the initial diagram, step N, the wake has just approached 
the cascade and is being cut apart by the blade. So far the 
only noticeable effect of activity of the vortices 15 the reduced 
velocity of their motion, compared to the potential flow 
velocity (compare the positions of the active vortices and the 
corresponding passive dots). The reduced velocity of motion 
of the von Karman vortex street in the flow, resulting from an 
additional velocity component induced by the street on itself, 
was first proved, both theoretically and experimentally, by von 
Karman [32], and its presence here testifies to the correctness 
of wake motion modelling. The next diagrams, steps N + 0.2T 
and N + 0.4T, show early stages of wake deformation. Negative 
vortices composing the central part of the wake segment are 
accelerated, while those situated close to its end are decelerated. 
An opposite tendency is observed in the motion of the positive 
vortices. As a result, the outer vortices of the two rows come 
closer to each other while the central ones become more distant. 
All this leads to noticeable segregation of vortices representing 
different rotation signs, with their further motion in two nearly 
isolated groups, visible in the final diagrams, steps N + 0.6T 
and N + 0.8T. Of these two groups, only the second group of 
vortices moves approximately with the velocity predicted by 
FlowER for the entire wake, while the first group moves faster, 
contributing to the growing rate of deformation of the inner 
structure of the wake. 

Beside the basic variant,  — 0, the study included three 
more variants of wake vortex distribution, which were 
denoted by phase shifts, Фф, equal to 0.25, 0.5, and 0.75. The 
above described segregation of a one-sign vortices and their 
further motion in nearly isolated groups were recorded in all 
cases. The most spectacular differences were observed in the 
motion of the vortices located close to the cut-off area. These 
vortices took visibly different trajectories, which in extreme 
cases led to the separation of isolated vortices from the wake. 
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Fig. 9. Stator wake deformation in the rotor passage: Г* = 0.005; o = 0.0 


Characteristic wake configurations with isolated vortices are 
shown in Fig. 10. 


Fig. 10. Extreme cases of wake deformation. Step N+T: separation of 
isolated positive vortex (left) and negative vortex (right) 


The above result prove and clearly illustrate the existence 
of mechanisms in the flow which make the same stator wake 
behave in a different way in the rotor passage only because of 
a different phase shift in the initial vortex distribution. In real 
wakes this distribution may be generally assumed a stochastic 
parameter, as for the time being, there is no evidence which 
would link the timing of the generation of stator wake 
vortices with the relative stator-rotor location. That means that 
consecutive segments of the same wake shed from a stator blade 
may behave in a different way in corresponding rotor passages. 
Or, in other words, consecutive segments of identical wakes 
shed from consecutive identical stator blades may behave 
in a different way in one and the same rotor passage. 


NUMERICAL AND EXPERIMENTAL 
VALIDATION 


The best way to validate the correctness of the above 
presented VDT results would be to compare them with the 
experiment, preferably performed in the same geometrical 
and flow conditions. Unfortunately, this way of validation is 
not possible in the present case, as the examined stage was 
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Fig. 11. Velocity fluctuations calculated for the passing vortex wake 
at selected points inside the rotor passage: dashed lines — individual 
realisations, continuous line — averaged realisation; ps, ss — rotor blade 
pressure and suction side, respectively 
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part of a real turbine in operation in a Polish power plant. 
Instead, an indirect method was applied which made use of 
the data recorded in other, but qualitatively similar conditions. 
It is noteworthy that despite a relatively large number of 
both numerical and experimental studies on W/R interaction 
which were recently published, finding appropriate data in 
which the wake development analysis goes as deep as to the 
vortex structure ofthe wake was extremely difficult. The most 
comprehensive and most valuable analysis from the point of 
view of the present validation was published by Kost et al 
[33]. Its authors presented the results obtained in a single- 
stage model turbine in the form of time-histories of velocity 
fluctuations recorded at points located inside the rotor passage 
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at selected dimensionless distances x/c from the rotor blade 
leading edge and at a distance equal to 1/16 of the rotor pitch 
from the blade surface, slightly off the rotor blade boundary 
layer. The experimental time-histories of velocity fluctuations 
presented by Koch et al. were averaged over 20 realisations. 
To provide better conditions for comparison, the VDT curves 
were also averaged, this time over 4 available realisations ofthe 
phase shift, ф = 0.00, 0.25, 0.50 and 0.75, as shown in Fig. 11 
in which the dashed lines represent individual realisations while 
the continuous lines show the result of averaging. 

Figure 12 compares the averaged VDT curves from Fig. 11 
with the time-histories presented by Kost et al. Due to different 
steady-flow velocity levels, the curves shown in the diagrams 
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Fig. 12. Velocity fluctuations inside the rotor passage: numerical results vs. experiment 
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represent velocity increments, dv, from an average level and 
not the absolute values. In each diagram the continuous line 
shows the averaged VDT time-history, while the two remaining 
curves represent the results presented by Kost et. al. The curve 
with thickened dots shows the experimental result and the 
dashed line - RANS calculations performed on 2D grid (only 
one realisation, without averaging). 

When comparing the experimental and VDT curves in 
successive diagrams we can notice high similarity of their 
patterns, especially at points located close to the pressure side 
of the rotor blade. Larger differences are observed at points 
situated close to the suction side, especially in diagrams for 
x/c = 0.1 and 0.7. Relatively high differences in the intensity of 
velocity fluctuations observed at point x/c — 0.1 may result from 
insufficient accuracy with which the vorticity distributed in the 
linear wake is modelled by a sequence of point vortices. When 
the wake deforms, these differences decrease significantly and 
only at point x/c = 0.7 remarkable differences can be observed 
again. This time they are of both quantitative scale referring to 
the amplitude of fluctuations, and qualitative nature observed in 
different shapes of the curves being the object of comparison. 
It 1s difficult to indicate clearly possible reasons of those 
differences, which would originate, for instance, from different 
exit flow conditions, which especially concerns flow angles 
calculated and recorded in the experiment. This suggestion is 
made more reliable by the fact that the results ofthe 2D RANS 
calculations reveal similar differences to those observed in the 
VDT time-histories. 


EFFECTS OF WAKE DEFORMATION IN 
VDT AND RANS FORMULATION 


The analysis described in the previous sections delivers 
a relatively clear and consistent pattern of the W/R interaction 
as an unsteady process controlled by the initial distribution of 
stator wake vortices at the time when the wake is cut apart by the 
rotor blades. Different distributions of these vortices provoke 
differences in their further behaviour and different effects of 
their interaction with each other and rotor blades. 

When confronting this interaction pattern with the RANS 
procedure of calculations, a basic question can be formulated 
whether and to which extent the RANS calculations, 
generally having the form of a single realisation, can be 
considered as representative for the general W/R interaction 
process. This question becomes even more important when 
we keep in mind that the RANS calculations, even performed 
in unsteady mode, provide rather limited opportunities for 
re-constructing the vortex structure of the stator wake inside 
the rotor cascade. Except special cases of integer stator/rotor 
blade number ratios, this reconstruction would require using the 
calculation region which covers all passages in stator and rotor 
rows. When we take into account the grid resolution required 
for proper capturing of the vortex interaction phenomena, 
meeting this conditions goes well beyond the calculating 
potential of even biggest computers in the world. In unsteady 
RANS calculations the above problem is usually solved using 
a special periodicity condition, bearing the name of phase 
lag, which allows the entire calculation region to be reduced 
to one stator and one rotor passage. Between the stator and 
rotor cascade a sliding plane is placed, on which the phase lag 
condition is met. When the wake is treated in a classical way, 
as a uniform structure shown in Fig. 4, the phase lag condition 
allows instantaneous wake parameters to be effectively passed 
between the stator and rotor calculation domains. Unfortunately, 
the vortex structure of the wake 15 lost in this process due to 
a conflict, in general, between the phase shift of the wake 


vortices and the phase lag defining the instantaneous position 
ofthe rotor with respect to the stator blades. As a consequence, 
the wake approaching the rotor calculation area is an averaged 
effect of various realisations representing different vortex phase 
shifts. In this situation the above formulated question can be 
rephrased to whether the wake *numerically averaged" at 
the inlet to the rotor calculation domain will behave in the 
same way, i.e. generate the same time-averaged effects as 
the real wake in multiple realisations. 


Fig. 13. Locations of recording points inside the rotor passage 


In order to find an answer to this question, comparison 
was made between velocity fluctuations generated by RANS 
and VDT codes at ten points situated close to the rotor blade 
trailing edge, Fig. 13. The dimensionless distances x/cofthese 
points from the rotor blade leading edge were equal to 0.1; 
0.3; 0.5; 0.7 and 0.9, while the distance from the blade surface 
was equal to 1/16 of the rotor pitch. The RANS calculations 
were performed in two variants differing by grid resolution. 
The first variant used the grid with 337 792 nodes (6032 
nodes in the circumferential plane y0z). This grid resolution 
can be considered similar to those used in turbine stage flow 
calculations performed for design purposes. The parameters of 
the second grid equalled to 698 880 nodes in the rotor passage 
(13 440 in the circumferential plane y0z) and were selected 
according to recommendation given by Arnone et al [18]. 

The results ofthe comparison are shown in Fig. 14. Dashed 
curves represent coarse-grid RANS calculations (variant 1) 
while those with large dots were obtained from the fine-grid 
calculations (variant 2). The third curve represents the VDT 
calculations used as the reference. 

As was shown in Fig. 9, the stator wake inside the rotor 
passage tends to take a trajectory close to the suction side of the 
rotor blade. That is why the velocity fluctuations recorded at 
the pressure side, points 2, 4, 6, 8 and 10, i.e. at a distance from 
the wake trajectories, are more regular and very similar to each 
other. What is also noticeable, the RANS results obtained for the 
fine grid are closer to the VDT curves than those obtained for 
the coarse grid. The observed differences between RANS and 
VDT curves seem to be accidental and do not reveal any regular 
tendency. Average velocities at the same points are slightly higher 
in RANS than in VDT calculations. 

In the vicinity of the suction side of the rotor blade, points 
1, 3, 5, 7 and 9, the differences in the shapes of the RANS and 
VDT curves are more pronounced. Paradoxically, the largest 
differences can be observed in the leading section located 
close to the rotor passage inlet. Downstream, the shapes of the 
corresponding curves become closer to each other. Like for the 
pressure side, the RANS results obtained on the fine grid are closer 
to the VDT results than those obtained on the coarse grid. 

Special attention should be paid to the results recorded at point 
9. This point 15 located near the rotor passage exit, close to the 
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Fig. 14. Velocity fluctuations inside the rotor passage: the effect of grid resolution 
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rotor blade suction side, and close to the hypothetical trajectories 
of vortices composing the stator wake. In these circumstances, 
all physical differences between the RANS and VDT approaches 
to modelling the W/R interaction, and all numerical errors that 
cumulate during the step-by-step calculations are expected to be 
the most remarkable at this point. However, the RANS curve 
obtained on the fine grid reveals a visibly similar pattern to that 
shown in the VDT calculations. Moreover, comparable locations 
and intensities of local extremes are also observed. 

The most remarkable difference between RANS and VDT 
results at point 9 is the average velocity level, visibly lower in the 
RANS calculations. The reason for this difference results from the 
fact that the boundary layer is thicker in this passage area and the 
recording point, unlike all other points, is situated inside it. 
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Fig. 15. Definition of velocity fluctuation amplitude 


The velocity fluctuation curves observed for coarse-grid 
RANS calculations suggest that in this case any quantitative 
assessment of the global effects ofthe W/R interaction on stage 
performance (forces, losses, efficiency) can be burdened with 
rather substantial error. 
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Fig. 16. Velocity fluctuation amplitudes recorded along the rotor passage 


To assess the correctness of modelling the rate of stator 
wake dissipation in the 3D RANS model, velocity fluctuation 
amplitudes were calculated at the rotor passage points from Fig. 13 
in the way illustrated in Fig. 15. The results of these calculations, 
shown in Fig. 16, reveal that increasing the grid resolution in 
RANS calculations leads to considerable reduction of wake 
dissipation. This tendency is illustrated by higher consistency of 
the amplitude change pattern obtained from RANS calculations 
with that representing the VDT solution, in which the wake 
dissipation is not directly taken into account. 

Good agreement between the VDT and RANS (fine grid) 
results shown in Figs. 14 and 16 is highly valued as a definitely 
positive answer to the question asked at the beginning of this 
chapter. The RANS calculations, performed on a sufficiently 
fine grid and making use of a “numerically averaged” stator 
wake, correctly predict time averaged effects of the real 
course of W/R interaction. 


CONCLUSIONS 


The article discusses the interaction between the stator 
wakes and the rotor cascade (the W/R interaction) in 
a turbine stage. Taking into account most recent result of both 
experimental and numerical studies of the inner structure of 
the wake, the stator wakes were modelled as segments of von 
Karman vortex streets. To achieve sufficiently high accuracy of 
unsteady calculations presenting consecutive stages ofthe W/R 
interaction, the rotor cascade was represented by an infinite row 
of blade profiles, the flow through which was calculated using 
the potential flow with vortex singularities and the conformal 
mapping method. Geometrical parameters of the stator wake 
and inlet conditions for the rotor flow were reconstructed from 
the preliminary 3D RANS calculations done for a stage, the 
geometry and flow conditions of which were taken from a real 
high pressure turbine in operation in a Polish power plant. The 
study has demonstrated that the deformation of the stator wake 
on its way through the rotor cascade leads to the segregation 
of two groups of one-sign vortices which then move down 
the passage along trajectories situated, as a rule, close to the 
rotor blade suction side. Different phase shifts in the initial 
distribution of vortices in the stator wake result in noticeable 
differences in the wake deformation course. 

The course of deformation predicted by the VDT analysis 
was compared to that observed experimentally and analysed 
numerically by other researchers in a turbine stage. The 
comparison revealed good qualitative agreement both with 
respect to particular stages of the wake deformation and its 
effects, having the form of velocity fluctuations recorded at 
selected points near the pressure and suction sides of the rotor 
blade. 

The last section of the article analyses the potential of 
RANS codes for determining quantitative effects of real 
W/R interaction in a turbine stage. The assumed criterion 
of applicability was the level of agreement of the RANS 
results with the time-averaged VDT data. It was stated that 
despite certain differences in the wake structure model used 
in each of these two codes and the resultant different nature of 
wake behaviour in the rotor passage: deterministic in RANS 
calculations and stochastic in VDT calculations, the RANS 
results obtained on a grid with a sufficiently fine resolution 
correctly represent the time-averaged effects of the real W/R 
interaction. The process of averaging of the vortex structure 
of the wake on the sliding plane between the stator and 
rotor calculation areas can be considered sort of “numerical 
averaging” of different realisations of the real vortex W/T 
interaction. 
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ABSTRACT 


This paper presents a methodology to evaluate the technical state ofa Francis turbine by shaft rotor dynamic 
simulation. There are several rotor dynamic criteria that define the technical state of a turbo-machine. 
To feed the shaft rotor dynamic model this delivers the required information to accomplish the technical 
assessment. The numerical rotor dynamic model uses as input, the field forces obtained by the fluid-solid 
interaction analysis undertaken over the blades of the runner. 
The rotor dynamic numerical simulations allow to determinate the record-in-time of the displacements of any 
point along the shaft. This information is relevant for diagnosis tasks, because it is possible to decompose 
it spectrally and to estimate the severity of the vibrations. Comparing the results of the numerical model 
against those obtained from machines that operates under normal conditions, it is possible to determinate 
the technical state of the turbo-machine. This allows studying the stability of the turbine working on 
several operation ranges. 
A Francis turbine is a very complex machine that involves many physical phenomena of different nature. 
In this way, the hydraulic input forces needed by the rotor dynamic model should not be assumed but 
calculated directly from the fluid interaction over the turbine structure. 


Key words: FEM; Francis turbine; diagnosis; fluctuations; hydrodynamic 


INTRODUCTION 


The technical diagnosis is a set organized of methods and 
a half that allow the evaluation of the technical condition 
state (his reason, evolution and consequence) of systems of 
engineering. In most cases they are systems in operation, 
designed to fulfill a principal function. These appliances 
generate and transform information, which is used later to study 
his technical condition state. 

For the definition of the field of application of the technical 
diagnosis, are studied the processes of destruction of the objects 
during the whole period of his existence and certain situations, 
in which it is necessary to consider the actions (shares) of 
diagnosis. The term "period of existence ofthe machines" refers 
to the cycle that initiates in the instant that there are formulated 
the requirements of design that must fulfill the new object up 
to his respective liquidation (residues, felling and reutilization) 
(Zóltowski B., 2002). 

The current technology (skill) is a result of the models’ 
analysis, which they describe with certain degree of accuracy. 
The process, which takes as an aim the construction of the 
best model of operation (mathematical or empirical), is 
named a process of identification, and includes shaped, you 


try, estimation and check. The methodology of the studies 
of diagnosis includes the heuristic and mathematical models 
(Żółtowski B., 1996). 

“The identification in the diagnosis is a process of searching 
complex, in that, based on the knowledge of the structure and 
the functioning of the object, there are realized the analysis 
and the synthesis of the model of the object for the obtaining 
of the information of diagnosis”. The idea on the methodical 
construction of the dynamic model of the mechanical system is 
supported of the information of the modal experiment and the 
information of the measurements of the quantities of vibration 
(acceleration of the vibration, moment rotational, force and 
etc.) in conditions of exploitation development. The models of 
analysis of loads are constructed depending on his application 
in the procedures of dynamic optimization and diagnosis of 
exploitation development. (Zóltowski B., Cempel C., 2004). 

During the electric generation process, the turbines 
develop hydraulic phenomena that extend along all the 
components of the machine represented in vibrations, sound 
and instabilities. Recently, the evaluation of the machine state 
using modern diagnosis techniques has become important. 
These diagnosis techniques, expert systems in particular, are 
based on the knowledge of the technical system and on the 


POLISH MARITIME RESEARCH, No 2/2013 41 


use of artificial intelligence. To enlarge the understanding of 
machines for failure detection, the study of their dynamical 
behavior plays a significant role. This methodology is known 
as maintenance based on models and pretends to achieve the 
failure identification by establishing relations between the real 
object and a virtual model (Kicinski, 2006). 

The vibration analysis has been a important tool in the 
technical diagnosis. Instruments for measuring vibrations 
(accelerometers) are implemented in large machines and 
the parameters extracted from the acquired signals are 
used as effective indicators of the technical state. The rotor 
dynamical study of turbines, in this case, pretends to develop 
the knowledge to correctly identify the relations between the 
symptom and its causes. 

The objective of the present study is to develop a study 
method that allows calculating the displacements of any point 
along the shaft of the turbine using FEM simulations. This 
information is relevant for diagnosis tasks, because it is possible 
to decompose it spectrally and to estimate the severity of the 
vibrations. 


CASE OF STUDY 


This investigation considers a Francis turbine composed by 
20 stay and guide blades and 15 runner blades. It rotates at 900 
rpm and generates a maximum output of 10MW. For this case, 
only the Best Efficiency Point (BEP) is analyzed. The turbine 
head is 230 m and the discharge at this point is 4.8 m/s. 

The assembly of the turbine / generator can be observed 
in Fig. 1. The set is assembled on a rigid shaft, to which are 
attached: the runner, the flywheel, the generator core and poles, 
and the exciter. The shaft is supported by three hydrodynamic 
bearings. The first of them (nearest to the runner) has axial and 
radial effect; the other two have only radial effect. 


METHOD OF EXPERIMENTAL MEASURES 


A series of measurements using dynamic pressure sensors 
arranged along the guide vanes and draft tube were carried 
out. Also, a set of strain gages installed in the turbine shaft 
were employed with the purpose of measuring deformations 
to determine the torque and the radial forces. 
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Runner (Francis Turbine) 


There were installed sensors along different components of 
the turbine. Three of them were installed near the zone where 
the rotor-stator interaction occurs (between the guide vane and 
the runner). Fig. 2 shows an image ofthe approximate position 
of one sensor. 


Fig. 2. Position of the sensors 


The task ofthese sensors was to capture pressure fluctuations 
generated by the rotor-stator interactions. Fig.3 shows an image 
ofthe experimental assembly carried out. It can be appreciated 
the physical conditions of the sensor installation. 


Fig. 3. Photo of the experimental location of the sensor 


Another six dynamic pressure sensors were installed along 
the draft tube. The objective of the sensors was to obtain 


Coupling 3 (Rigid) 


Fig. 1. Turbine shaft with bearings and attached mechanical elements 
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information on the vortex rope generated at the exit of the 
runner. The sensors were distributed in two cross sections: i.e. 
A and B; in each section were installed 3 sensors (Fig. 4). 


SECCION B 


SECCION A 


Sensor 1 


Sensor 2 Sensor 3 


SECCION 


Fig. 4. Distribution of the sensors in the draft tube 


The measurements were recorded performing a simultaneous 
sampling during 3 seconds of every channel at 102 kHz, i.e. 
more than 40 rotor revolution and more than 6500 samples 
per channel. It allows to capture all the hydraulic phenomena 
present at low and high frequencies. 

For the measurement of the torque and the radial force 
a strain gage technique applied to the generation shaft was used. 
The use of this technique did not affect the turbine operation in 
any way. The strain gage signals were registered via wireless 
connection. 

Fig. 5 shows the position in the shaft where the strain 
gages (SG) were installed. For the measurement of the radial 


force it was use of four bi axial well doable SG meanwhile 
the measurement of the torque used one bi axial well doable 
SG. The strain gages were connected to a transmitter tied to 
the shaft. 


Bearing 


Runner flange 
Fig. 5. Strain gages position in the shaft 
Fig. 6 shows an image of the experimental installation of 


the SG. The distribution of the strain gages and its position in 
the shaft can be observed. 


Ји 


Fig. 6. Photo of the experimental installation of the gages in the shaft 


The information give by all this pressure sensors and the 
strain gages is used to validate the CFD numerical model. 


ROTODYNAMIC METHOD 


The goal is to develop a model that accurately predicts the 
dynamic behavior of the shaft of a Francis turbine by means 
of numerical methods. To achieve this goal, it is necessary to 
find an approximated solution to the equations of motion of 
the system. These equations are a set of non linear differential 
equations: 


Mx + р(х, i) + K(x,3)x = P(t) (1) 
where: 
M  — global matrix of inertia 
D  -— global matrix of damping 
К — global matrix of stiffness 
P(t) — generalized vector of external excitations 
X  — generalized vector of position 
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The model is implemented in the computational environment 
MESWIR®. This software was developed by the Institute of 
Fluid Flow Machinery (IFFM) from the Polish Academy of 
Science. It composed by a set of codes to determine static 
and dynamic parameters of rotors, including analysis of the 
hydrodynamic bearings. The software uses the finite element 
method to analyze the rotor shaft and the finite difference 
method to determine the characteristic of the hydrodynamic 
bearings. With the software it is possible to analyze lateral 
vibrations with linear and non linear models. 

The required entries of the model are: the conditions of 
support, the forces or external excitations, the mechanical and 
geometrical properties of the shaft. The shaft is supported in 
three hydrodynamic bearings but for the rotordynamic model 
they are represented by springs and dampers. The determination 
of the oil film stiffness and damping coefficients required 
additional computation of the fluid dynamics involved. The oil 
film stiffness and damping coefficients depend on the relative 
position and velocity of the journal in the bearing causing the 
non linearity in the equations of motion. However the results 
obtained under the assumption of linear stiffness and damping 
coefficients are acceptable. For the purpose of this work, the 
MESWIR® linear models were used. 

The bearings are also attached to a foundation. The 
foundation is replaced in the mechanical model by springs 
and dampers, whose coefficients were prior determined by 
experimental measurements and numerical simulations. The 
shaft is modeled coupled rigid discs with masses and inertias 
equivalent to those of the real objects. 


Geometry and physical properties 


The real model is simplified in a mechanical model to make 
possible its analysis (Fig. 7). The shaft is discredited into 33 
cylindrical Timoshenko beam elements with six degrees of 
freedom to use the finite element method. These elements are 
provided of mass, stiffness and damping, which are estimated 
from the physical properties of the material of the shaft, such 
as: density, Young modulus, Poisson coefficient and damping 
coefficients. In the Tab. | are listed the material properties. 


Unbalanced mass 1 


Hydraulic Forces 


Bearing 1 С] 


tunner 


Foundation 


Flywheel 


lac Bearing 2 


Tab. 1. Physical properties of the material of the shaft. 


Young i Damping 
coefficient | coefficient A 


Damping 
ys 


7800 5.09*1055 


Support parameters 


In the oil film in the hydrodynamic bearings appears 
a hydrodynamic pressure distribution induced by the rotation 
and originating the reaction force. This pressure distribution 
18 described by the Reynolds equation: 
= (2) 
ot 


ee 2). eia 
R? 204 00) Oz Oz R 60 


where: 

— pressure in the oil film 

— film thickness 

oil velocity in tangential direction 
— dynamic oil viscosity 

— radius of journal 


юз с = 
| 


Assuming that the flux is laminar and that the viscosity is 
constant along the oil film the pressure distribution could be 
found. The magnitude and orientation of the reaction force 
depend of the position and velocity of the gravity center of the 
rotor shaft. The stiffness and damping coefficients are computed 
based on the reaction force. The stiffness coefficients are the 
derivatives of the reaction force respect to the displacements. 


OW OW. 
kg = =, ‘ky = : 
Ox ду 
3 
oW, oW, " 
Ку = > Ко > 
ду дх 
where: 
Ку, ky», Ко, Кк = oil film stiffness coefficients 
Wo W, — components of the reaction force 


Unbalanced mass 2 


• 
4 Unbalanced mass 3 


Electromagnetic Forces 


Exciter 


ij 
Bearing 3 Le] 


Rotor of the generator 


Foundation Foundation 


Fig. 7. Numerical model 
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The damping coefficients are found as the derivatives of 
the reaction force respect to the velocities. 


e ow, -— oW, 
11 ГТ > 22 oy а 
OW. OW 


where: 
Су, Сур Сүз, Су — ОП film damping coefficients 


The Fig. 8 presents a representation ofthe oil film as springs 
and dampers. 


e 


Fig. 8. Stiffness and damping coefficients the oil film 


The shaft is affected by static and dynamic loads. The static 
loads are the masses of the elements attached to the shaft. The 
dynamic loads are originated by: residual unbalanced masses, 
the solid-fluid interaction in the runner and electromagnetic 
phenomena in the generator. 


Static and dynamic loads 


The static loads are determined by the geometry and 
physical properties of the discs assembled to the shaft. For 
the model, it is supposed that these discs are made of steel. 
The discs add mass and inertia to the system but they do not 
represent additional stiffness or damping. Other properties of 
the discs are list in Tab. 2. 


Tab. 2. Elements assembled to the shaft 


Outer 
diameter [m] 


Generator 


poles Exciter 


Flywheel 


Inner diameter 
[m] 
Thickness 
[m] 

Mass 
[kg] 
Inertia 
[kgm?] 

As it was mentioned the dynamic loads could be inertial, 
hydraulic or magnetic. The inertial forces are used to tune 
the model and the represent an output of the model. The 
hydraulic forces are estimated by CFD simulations. Regarding 
to the magnetic forces only unbalanced magnetic pull is 
considered. 

The unbalanced magnetic pull appears when the rotor 
is eccentrically positioned respect to the stator. This causes 


asymmetric distribution of the magnetic flux and then an 
unbalance force in the direction of the smallest air gap. This 
force depends in a non linear way of the relative eccentricity, 
which is defined as: 


£ - u/AR (5) 
where: 
=  — relative eccentricity 
u — distance between the rotor center line and the stator 


center line 
AR - difference between the stator inner radius (R,) and the 
rotor radius (К) 


The unbalance magnetic pull has a constant and a variable 
part. The variable force oscillates with a frequency equal to the 
supply frequency. If the number of poles increases, the variable 
part decreases. The mean value for the unbalance magnetic pull 
could be found as (Gustavsson, 2005): 


цоЅ26 Ьл е 


=== 6 
2p2AR2 mmc (6) 


where: 

Ho = magnetic permeability 
5, = lineal density of current 
R, — stator inner radius 

h – rotor length 

e — relative eccentricity 

p - nNumber of poles 


By means ofa static analysis the deflection of the shaft in 
the generator core section 15 found. With this value and the 
technical information of the machine the unbalance magnetic 
pull is computed and is found that it has a value of 4966 N. 


RESULTS 


The experimental data used in rotor dynamic calculations is 
obtained from displacement sensors installed in the bearings. 
In each bearing are two sensors forming a rectangular angle 
between them and an angle of 45 degrees respect the horizontal 
plane as is shown in Fig. 9. Then, the vibration is measure 
in a coordinate system that is 45 degrees rotated respect the 
general coordinate system and some transformations should 
be made to compare the results. 


Fig. 9. Place of the sensors on the rotor 


The most relevant component in the spectrum is at 15 Hz, 
which is the rotational frequency of the machine. The other force 
components are very low and are not considered in this analysis. 
In the model the force should be introduced as an unbalanced 
mass at specific unbalance radius. The transformation is made 
using the definition of centrifugal acceleration. 


POLISH MARITIME RESEARCH, No 2/2013 45 


— —- Horizontal Force 


700 * = * Vertical Force 


0 - = а аса 


0 50 100 150 200 250 300 350 400 
Frequency [Hz] 
Fig. 10. Spectral analysis of the radial force 


Tab. 3. Spectral components of the radial force. 


Frequency Horizontal force Vertical force 


[Hz] IN] IN] 


Horizontal displacement 


т 
e 


> 
© 


Displacement [um] 


-30 


0.02 0.03 0.04 0.05 0.06 0.07 0.08 
Типе [5] 


Displacement | um] 


0.02 0.03 0.04 0.05 0.06 0.07 0.08 


Time [s] 


Fig. 11. Numerical and experimental displacement signals for the bearing 1 
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The analyses are made with the modules of MESWIR® for 
studying lateral vibrations with linear models. Hence, only 
excitations at first harmonic are taken into account and the 
numerical signals are sine functions with the rotational frequency. 
In Fig. 11, 12 and 13 the experimental and numerical results for 
the displacement vibration are compared for each bearing. 

Horizontal displacement 


25 


Displacement [um] 
tn 


0.03 0.04 0.05 0.06 
Time [s] 


Displacement [um] 


-40) 
0.03 0.04 0.05 0.06 0.07 0.08 0.09 
Time [s] 
Fig. 12. Numerical and experimental displacement signals for the bearing 2 
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-10 


-30 


0.05 0.06 0.07 0.08 0.09 0.1 


Time [s] 
Fig. 13. Numerical and experimental displacement signals for the bearing 3 


As it was mentioned the unbalanced forces are tuned until 
the numerical vibration amplitudes adjust to those measured 
in the real technical object. The final required values represent 
a result from the rotor dynamic model. 


Tab. 4. Resultant unbalance masses 


Unbalance | Unbalance 
d mass e radius 


[kg] [m] 


Phase 


Position [degrees] 


Runner 


Flywheel 


Generator poles 


Runner 
(hydraulic force) 


Runner 
(hydraulic force) 


FINAL CONSIDERATIONS 


The radial forces obtained via numerical methods were 
used to feed the FEM rotor dynamic software. These forces 
have a consistent predominant spectrum component at the 
frequency of rotation which was used in the linear rotor 
dynamic model. However this resultant radial component 
has a low magnitude, because for the simulations the runner 
is completely symmetrical and the unbalance in the pressure 
distribution is little. If some imperfections of the runner are 
included in the simulations, the magnitude of force would 
probably become more significant. 

It was possible to develop a rotor dynamic model whose 
results are in agreement with the experimental data. In this 
case, the unbalance masses required are relative large. This 
could be caused because the rotor is very rigid and important 
amounts of residual unbalanced could be bear with acceptable 
vibration amplitudes. 

The amplitudes of vibration for the displacement obtained 
by the numerical computations and the experimental measures 
are acceptable according to the standard ISO 7919-5. This 
standard establishes a maximum limit of 165шт. From the 
point of view of vibration the technical state of the machine 
is acceptable. 
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Design analysis of ORC micro-turbines making 
use of thermal energy of oceans 


Marian Piwowarski, Ph.D. 
Gdansk University of Technology, Poland 


ABSTRACT 


The article presents the results of the analysis of energy conversion cycles making use 
of thermal energy of oceans. The objects of analysis were two cases of closed Organic 
Rankine Cycle (ORC) power plants, which were: the cycle in which the vapour of the 
working medium was produced by warm oceanic water in the circum-equatorial zone, and 
the so-called "arctic" cycle in which this vapour was produced by non-frozen water in the 
circumpolar zone. Between ten and twenty low-boiling media were examined for which 
operating parameters were optimised to obtain the highest cycle efficiency. A preliminary 


design of an ORC turbine which was obtained by optimising basic design parameters is included. It has 

been proved that realisation of the Ocean Thermal Energy Conversion (OTEC) cycle is possible both in the 

warm and permanently frozen regions. The results of the calculations have also revealed that the efficiency 

of the OTEC cycle is higher in the circumpolar zone. Selecting a low-boiling medium and designing a highly 
efficient turbine operating in both abovementioned regimes is technically realisable. 


Key words: Ocean Thermal Energy Conversion (OTEC) cycles; Organic Rankine Cycle (ORC); 
low-boiling media; steam turbines; turbine design 


INTRODUCTION 


Seas and oceans are the world largest reservoirs in which 
thermal energy of sunbeams is accumulated. The sunlight falling 
on the ocean surface is intensively absorbed by the water which 
collects all its energy in the shallow surface layer. Unlike the 
solar energy, the intensity and concentration of which is subject 
to steady cycles and accidental disturbances, the energy of seas 
and oceans reveals stable intensity, both in 24 hr and annual 
time scale. This is an important advantage as it leaves aside 
the problem of energy storage and its costs. In tropical regions 
of oceans, approximately between 15? north and 15? south 
latitude, the heat collected from sunbeams heats the water in 
the shallow surface layer to the practically constant temperature 
nearing 28 °С. The annual average temperature in the entire 
region ranges approximately between 27 °С and 29 °С. Below 
the warm surface layer the water becomes colder, and at the 
depth between 800 and 1000 m the water temperature drops 
down to about 4 °С. Depending on the region, temperature 
changes can take different courses and for instance, at the 
depth of 200 m the water temperature can be already as low 
as 7 ?C. The temperature difference is kept constant all year 
round, with small deviations of a few Celsius degrees resulting 
from seasonal weather and day/night changes. Utilisation of 
this difference bases on the application of a working medium 
which evaporates at the surface temperature and condensates 
when cooled by the water taken from the much deeper level. 
This type of energy conversion is frequently referred to as 
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Ocean Thermal Energy Conversion (OTEC). A key issue 
in case of localisation of a number of OTEC units close to 
each other is the distance between them. Due to outlet water 
disposal, this distance should not be excessively small. This 
aspect is of high importance as the surface waters are cooled 
each time the outlet water is disposed from an open OTEC 
installation. In the system of combined units the excessive 
amount of water disposed by the system can considerably 
decrease the temperature of the surrounding regions, which 
in turn can affect the efficiency of all power plants situated 
nearby. Initial concepts on how to make use of the energy 
collected in seas were presented by J. d'Arsonval (France) in 
1881. The first experimental object ofthis type was built by his 
pupil, G. Claude, in Matanzas near Cuba in 1928. The power 
output of the installation was 22 kW and it worked only for 
2 weeks, after which was partially destroyed by typhoons. In 
late fifties of the last century French researchers began again 
to study the subject, building an installation in the Ivory Coast 
region. Further development in this field was initiated in 1960 
by Andersens, two Americans, father and son, who established 
a private company exploring thermal energy from seas using 
OTEC technology. А few years later a 50 kW semi-technical 
system was built in the USA, which was followed by a fully 
technical one of 1 MW power output, built in 1980 and installed 
on a reconstructed ship. Cold water amounting to 450 tonnes 
was delivered using polyethylene pipes having the diameter of 
1.5 m and length of 700 m. In 1981, on Nauru island in Japan 
a stationary 100 kW installation was built, to which 1410 t/h 


Tab. I. List of demonstrative OTEC power plants [8] 


Location 
Cuba 


- 


1930 


Producer 
Claude, France 


Power [kW] 
Open 


1979 Hawaii 


Mini OTEC, USA 


Closed 


1980 Hawaii 


OTEC-1, USA 


Closed 


1982 Nauru 


Toshiba& TEPC, Japan 


Closed 


1984 Saga 


Saga University, Japan 


Open 


1992 Hawaii 


Nelha, USA 


Closed 


1992-1998| Hawaii 


о | 5 | 2 | ел] 5 | 6 | 


Nelha, USA 


Closed 


1995 Saga Saga University, Japan 9 Closed 


© 


2000 Tuticorin 


of cold water having the temperature of 8 °С was delivered 
from the sea bottom using a 1093 m long pipe. In this area the 
temperature of surface waters nears 30 °С. About 75 % of the 
produced energy was used for own purposes of the installation 
and the remaining part was transmitted to the electric network. 
Tab. 1 collects demonstrative OTEC power plants which were 
the objects of examination. 

Plans are to build a power plant fed by the Gulf Stream 
which flows near the USA. The temperature difference between 
the Gulf Stream and the surrounding waters secures efficient 
recovery of thermal energy at remarkable cost reduction 
resulting from a short distance between the lower and upper 
heat sources, and the disposal of the used water to the Gulf 
Stream which will transport it away. At the beginning of the 
21-th century, numerous countries, mainly Australia, Germany, 
France, the United Kingdom, the Netherlands, Italy, India, Japan, 
Canada, and the USA, began to value thermal energy collected 
from the oceans to produce electric energy and intensified their 
activities oriented on development of technologies leading to 
most effective use of this energy [2+7, 9, 13, 15, 18, 20-22]. 
Technical solutions of these installations base on the Open Cycle 
- Ocean Thermal Energy Conversion (OC - OTEC) cycle and 
the Closed Cycle — Ocean Thermal Energy Conversion (CC 
— OTEC) cycle. Open systems make use of the water pumped 
directly from the ocean. The water evaporates when flowing 
through the low-pressure chamber in which the pressure is 
decreased below the condensation point. This way salt is 
separated from the water. Thermal energy needed for steam 


NIOT, Indie 


1000 Closed 


production in the evaporator is released by the surface water 
which is the upper heat source. After leaving the evaporator, 
the produced steam expands in the turbine and then flows to 
the condenser cooled with the water taken from deeper levels 
of the ocean. Oxygen, carbon dioxide, and nitrogen, which 
have been dissolved in the oceanic water and are not subject to 
condensation in the condenser, are released to the environment 
using the vacuum pump which increases the pressure to the level 
required for purifying the installation from these gases (Fig. 1). 
A characteristic feature of open installations 15 relatively high 
speed of medium evaporation, which is directly related with 
the processes of mass and heat transfer. The steam produced 
in the heat exchanger reveals low, or zero salt concentration, 
along with relatively low toxicity. Further possible utilisation of 
the medium used in the system can be an additional attractive 
feature of the discussed cycle in places where clean water 1s 
hardly accessible and highly required. 

In case of closed systems, the principle of their operation 
is the same as in the Rankine cycle (Fig. 2). In this concept the 
surface water is also used as the upper heat source. The energy 
from outside 15 passed to the medium in the heat exchanger 
to make it evaporate. Due to smaller specific volume of the 
working medium, the dimensions of the devices are smaller 
than in open cycle installations. Proposals can be found in the 
literature to combine OTEC systems with solar collectors, thus 
creating so called Solar and Ocean Thermal Energy Conversion 
(SOTEC) systems. Such solutions increase the efficiency ofthe 
entire cycle [6], [22]. 


Steam 


Steam 
generator 


Generator 


Condenser 


Cold water 


Fig. 1. Scheme of open OTEC cycle 
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Fig. 2. Scheme of closed OTEC cycle 


In each cycle making use of thermal energy of seas 
and oceans, particular attention 1s paid to leak tightness of 
installation components. Due to huge mass flow rates and high 
requirements for warm and cold water, remarkable part of the 
produced electric energy is used for covering own needs of the 
installation, mainly connected with pump feeding. Elimination 
of leakages to the environment improves the overall efficiency. 
Both organic and non-organic media can be used in CC - OTEC 
cycles. Low-boiling media ideally meet conditions for liquids 
used in those units, and their biggest advantages are low 
boiling pressure and temperature which is frequently lower 
than 0°. However, their use is sometimes connected with certain 
inconveniences, as they reveal, to a smaller or greater extent, 
flammability and/or toxicity. Their possible direct contact with 
the environment is highly inadvisable. Advanced protection 
concerning safety of work, i.e. leakage tightness of the closed 
cycle, is in this context fully justified. 


LOW-BOILING MEDIA 


Working media should meet various requirements which 
are not connected with their thermal characteristics. First of all, 
they should not be toxic and should be environment friendly. 
They should be cheap and easily accessible. They cannot react 
with flow system components. А properly selected low-boiling 
medium should provide opportunities for system operation at 
highest possible efficiency, and for maximal possible utilisation 
ofthe available heat source. Other very important criteria which 
are to be met by the working medium and are to be taken into 
account in its selection include: low freezing temperature, high 
critical temperature, corrosion passiveness, lack of toxicity, 
incombustibility, low price, etc. Moreover, it is important 
for the working medium not to produce any threat for human 
beings and the environment. The effect of a given medium on 
the environment is measured using ecological coefficients, 
mainly ODP (Ozone Depletion Potential) and GWP (Global 
Warming Potential). 


ASSESSMENT OF RESULTS OF 
CALCULATIONS OF CIRCUM- 
EQUATORIAL OTEC CYCLE 


Within the framework of the present analysis, calculations 
were performed for the closed OTEC cycle working in 
conditions characteristic for circum-equatorial regions (see 
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Fig. 2). The temperatures of the media assumed for the 
calculations are shown in the figure. The next assumption was 
the mass flow rate of 100 kg/s used for vapour production in 
the evaporator. For the closed OTEC cycle power plant several 
tens of media were analysed, after which more than ten media 
securing high cycle efficiency and technical realisability were 
selected [10, 11, 19]. The below presented results refer to the 
following selected media: Water, Acetone, Methanol, Isobutane, 
Dimethyl ether, R114, R227ea, R245fa, Trifluoroiodomethane, 
R124, R1234ze, Ammonia. For each medium the pressure at 
turbine inlet was optimised in such a way that maximum power 
was reached, and the heat exchanger design was technically 
realisable. 

From the warm water mass flow rate equal to 100 kg/s we 
can produce, approximately, from 48 to 50 kW of turbine shaft 
power (see Fig. 3), which corresponds to cycle efficiency not 
exceeding 2.9 % (see Tab. 2). The power outputs obtained in the 
analysis were close to each other, therefore medium selection 
should be based on other parameters which would secure 
designing a cheap, reliable, and highly efficient installation. 
The experience gained in Japanese experiments suggests that 
the costs of the working medium, turbine, pumps and heat 
exchangers amount to as much as about 20 %. In this situation 
attempts are to be made for the turbine cycle units to be properly 
optimised. Main attention was focused here on optimisation 
of the design parameters of the turbine working with selected 
low-boiling media [14]. 

The analysis included optimisation of design parameters 
for a single-stage axial turbine with respect to the generated 
power output. The optimised parameters were the following: 
velocity coefficient, reaction, rotational speed, and supply 
arc dimension. An attempt was also made to keep the ratio 
of the effective diameter to the blade length larger than 4. 
Another parameter which was also assessed and taken into 
account in eth optimisation was the Mach number. Tab. 1 
presents a collection of optimal parameters for the single- 
stage axial turbine. In steam turbines typical values of the 
velocity coefficient depend on the type of the turbine stage 
and range from 0.4 to 0.5 for impulse stages, and even as much 
as from about 0.6 to 0.7 for reaction stages. For the analysed 
media the selected velocity coefficient not exceeded 0.6 (see 
Fig. 4), which is a typical range for systems of this type. The 
stage reaction was selected in such a way that the peripheral 
efficiency was the highest possible, keeping at the same time 
the hub reaction not larger than 0.05. Typical reaction used in 
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Fig. 3. Turbine shaft power for the examined media 


Fig. 4. Single-stage turbine velocity coefficient for the examined media 
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Fig. 5. Single-stage turbine reaction for the examined media 
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turbine stages does not exceed 0.7. During the present analysis 
the stage reaction did not exceed 0.5 (see Fig. 5), i.e. it also took 
typical values [12]. The rotational speed obtained in the analysis 
depended on the type of medium and was equal to 3000 and 
6000 rpm for most media, while for Trifluoroiodomethane and 
Ammonia it was equal to 15000 and 16000 rpm, respectively 
(see Fig. 6). All these values are acceptable for technical 
realisation. In most cases full supply arc was reached, only 
for Trifluoroiodomethane it did not exceed 5 % (see Table 2). 
Reduction of the supply arc dimension increases remarkably 
the ventilation loss and decreases the internal stage work. This 
result, combined with slightly lower efficiency eliminates the 
above medium as a working medium. 

The turbine was designed in such a way as to obtain the 
length of the blades not smaller than 10 mm and the D/L ratio 
not lower than 4 (see Fig. 7 and Fig. 8). This condition was met 


for the reported media. The selected design parameters should 
have secured not exceeding 1 by the Mach number at stator 
and rotor exits. This condition was also met, and in each case 
subsonic stages were obtained (see Fig. 9 and Fig. 10). The 
performed analysis has proved that a highly efficient single- 
stage axial turbine can be designed for most of the examined 
media (excluding Trifluoroiodomethane). However, it seems 
that very low pressure (high negative pressure) eliminates 
water, acetone and methanol. In this situation, the media 
which looks most favourable from the technical point of view 
are R114, R245fa and Ammonia. Here, other criteria are to be 
taken into account, including toxicity of the medium (safety in 
case of leakage) and its effect on the environment. A detailed 
and realistic economic analysis is also to be performed. Final 
selection of a relevant medium will be a compromise between 
the advantages and disadvantages of this medium. 
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Fig. 6. Single-stage turbine rotational speed for the examined media 


Tab. 2. Comparison of calculated results for closed OTEC power plant with turbine design parameters; where: N jem orc-effective power of the cycle; 
Праг Cycle efficiency; Nearno Carnot cycle efficiency; рутефит vapour pressure at turbine inlet; Moa, medium vapour mass flow rate; H,-isentropic enthalpy 
drop in the turbine; v-velocity coefficient; p-stage reaction; Ma,,- Mach number at stage stator exit; ¢-supply arc dimension; n-rotational speed; D -mean 
stage diameter; l,-stator blade length; D/I-mean diameter-to-blade length ratio; Ma, ,- Mach number at stage rotor exit 


Parameter Nearnot H, 
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Medium/ 
Unit 
Water 


[kW] [kJ/kg] 
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Methanol 


Isobutane 


Dimethyl 
ether 


R114 
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R124 
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Fig. 9. Absolute Mach number behind single-stage turbine stator for the examined media 
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Ma, [-] 


Fig. 10. Relative Mach number behind single-stage turbine rotor 


RESULTS OF CALCULATIONS FOR 
ARCTIC OTEC CYCLE AND THEIR 
EVALUATION 


For this cycle a system was proposed which makes use of 
temperature difference between the frosty arctic air and the non- 
frozen oceanic water, which eliminates the use of long pipes 
for its collection. In the cycle shown in Fig. 11 a condenser 
cooled with frosty air was applied. Other media revealing low 
freezing temperature, aqueous solution of calcium chloride for 
instance, can also be applied [1]. 

The temperatures of the media assumed in arctic OTEC 
cycle calculations are shown in Fig. 11. The mass flow rate of 
the non-frozen water which was used for generating medium 
vapour in the vapour generator was assumed equal to 100 kg/s. 
In this case also several tens of media were analysed and, like 
previously, between ten and twenty media were selected as those 
which secure high cycle efficiency and technical realisability. 
Below are presented results of the analysis referring to the 
selected media, which were: Trifluoroiodomethane, Dimethyl 
ether, CO,, R23, R13, R218, Sulfur hexafluoride, R116, 
R245fa, R143a, Ethane, R124, R1234ze and Ammonia. Like 
in previous case, for each medium the turbine inlet pressure 


was optimised in such a way that maximum shaft power 
was obtained for this pressure and the design of the heat 
exchanger was technically realisable. Approximately, from 30 
to 80 kW of the turbine shaft power can be produced from the 
assumed mass flow rate of the warm water (see Fig. 12). In 
the examined arctic variant the obtained power levels differed 
much between each other. The largest power was produced 
using Trifluoroiodomethane, Dimethyl] ether, and Ammonia as 
working medium. As for the cycle efficiency, the above list of 
three media was complemented by R218, R245fa, R124 and 
R1234ze, for which the recorded efficiency exceeded 6 % (see 
Tab. 3). In the next step, optimisation of the design parameters 
for a single-stage axial microturbine was carried out. The list of 
optimised parameters included: velocity coefficient, reaction, 
rotational speed, and supply arc dimension. An attempt was 
made to keep the ratio of the effective diameter to the rotor blade 
length larger than 4. A parameter which was also assessed and 
taken into account in the optimisation was the Mach number. 
Table 3 presents a collection of optimal parameters for the 
single-stage axial turbine. 

The velocity coefficient selected for the analysed media 
ranged up to 0.7 (see Fig. 13). The stage reaction did not exceed 
0.5 (see Fig. 14). During the analyses different rotational speeds 
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Generator 
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————— 


Main pump 
Fig. 11. Scheme of closed OTEC cycle in arctic version 


54 POLISH MARITIME RESEARCH, No 2/2013 
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Fig. 12. Turbine shaft power for the examined media 


were obtained for different media, they took values from 3000 
rpm, to even as much as 37000 rpm for Ammonia (see Fig. 15), 
which was still technically realisable. In most cases a full 
supply arc was obtained, excluding Trifluoroiodomethane, 
Dimethyl ether, Ammonia and Ethane, for which the supply 
arc dimension did not exceed 20 % (see Table 3). In the latter 
situation the ventilation loss increases remarkable and the 
internal stage work decreases. The turbine was designed in such 
a way as to obtain the blade length not smaller than 10 mm 
and the D/L ratio not smaller than 4 (see Fig. 16 and Fig. 17). 
This condition was successfully met for the presented media. 
The design parameters were selected assuming that the Mach 
number at turbine stator and rotor exits does not exceed 1. This 


condition was also successfully met, and each time a subsonic 
stage was obtained (see Fig. 18 and Fig. 19). The performed 
analysis has confirmed that designing a highly efficient single- 
stage axial turbine 15 possible for the majority of the analysed 
media, from which: Trifluoroiodomethane, Dimethyl ether, 
Ammonia, R218, R245fa, R124 and R1234ze seem to be 
most favourable from the thermodynamic and technical point 
of view. If we take into account toxicity of particular media 
and their effect on the environment (see Tab. 4), the most 
favourable medium is R245fa which is non-flammable and 
nontoxic. However, a detailed and realistic economic analysis 
Is to be performed to select finally an optimal medium for the 
system of this type. 


Tab. 3. Comparison of calculation results for arctic OTEC power plant with turbine design parameters; where: Nem orc-effective power of the cycle; 
Праг CYCLE efficiency; Nanon Carnot cycle efficiency; py medium vapour pressure at turbine inlet; m),,,,-medium vapour mass flow rate; H,-isentropic enthalpy 
drop in the turbine; v-velocity coefficient; p-stage reaction; Ma,,- Mach number at stage stator exit; e-supply arc dimension; n-rotational speed; D -mean 
stage diameter; l,-stator blade length; D/I-mean diameter-to-blade length ratio; Ma ,- Mach number at stage rotor exit 


Parameter Tlcarnot H, 


Nas onc 


Medium/Unit | [kW] [-] 


Trifluoroiodo- 
methane 


[kJ/kg] 


Dimethyl ether 


Sulphur 
hexafluoride 


R116 


R245fa 


R143a 


Ethane 


R124 
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Fig. 13. Single-stage turbine velocity coefficient for the examined media 
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14. Single-stage turbine reaction for the examined media 
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Fig. 15. Single-stage turbine rotational speed for the examined media 
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Fig. 16. Single-stage turbine stator blade length for the examined media 


Fig. 17. Single-stage turbine D/L ratio for the examined media 


Fig. 18. Absolute Mach number behind single-stage turbine stator for the examined media 
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Fig. 19. Relative Mach number behind single-stage turbine rotor for the examined media 
CONCLUSIONS 3. Breeze P. Power Generation Technologies, ISBN 0 7506 6313 


To sum up, OTEC power plants are the source of clean 
and renewable electric energy. Along with the energy they 
can produce fresh water, which is of high importance in 
the areas with its shortages. A properly designed system 
generates limited, or even zero amount of carbon dioxide. It is 
optimistically estimated that the resources of thermal energy 
are so high that they can meet all current power requirements 
of the human society. Its disadvantage still remains high 
production cost, which is remarkably higher than in case of 
energy production from coal. For the operation of this type of 
power plant to be possible, the temperature difference between 
the upper and lower heat source should reach about 20 °C, the 
least. Moreover, vast resources of cold water should be situated 
close to the land for large power plants to be erected in the 
vicinity of human agglomerations. Up to the present, certain 
problems can be observed with collecting financial support 
for such actions, as only small-scale systems have been tested 
so far. Another disadvantage is unfavourable effect of on the 
environment, as localisation of power plants (pipes with warm 
and cold water) close to the land unfavourably affects the off- 
shore marine ecosystem and reefs. The here reported analysis 
has proved that constructing a system which utilises thermal 
energy stored in oceans 15 possible. It has also turned out that 
higher potential 1s this field can be found in the circumpolar 
regions. For the assumed mass flow rate of the water used for 
producing the organic medium vapour, certain media were 
proposed which secured installation efficiency exceeding 6 %. 
This result does not seem impressive, but the source of this 
energy is practically inexhaustible. 

For a group of low-boiling media a preliminary design 
of the axial turbine was presented. The realistic ranges of 
the selected design parameters have proved that constructing 
a highly efficient turbine making use of low-boiling medium 
Is possible. 
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ABSTRACT 


The aim of the present paper is to develop a computational fluid dynamics (CFD) analysis to study the 
combustion process in a four-stroke marine diesel engine, the Wärtsilä 6L 46. The motivation comes from 
the importance of emissions from marine engines in the global emissions, particularly for nitrogen oxides 
(NO) and sulfur oxides (SO). The pressure and temperature fields were obtained, as well as the exhaust 
gas composition. In order to validate this work, the numerical results were satisfactory compared with 
experimental ones, which indicates that this model is accurate enough to reproduce the fluid pattern inside 
the cylinder during the combustion process. Accordingly, the aim of future works is to use this numerical 
procedure to optimize the performance and reduce the emissions of the new marine engine designs. 


Key words: Marine engine; emissions; combustion; CFD (Computational Fluid Dynamics) 


INTRODUCTION 


Nowadays, the world's ships are mainly powered by diesel 
engines. Marine diesel engines are very efficient, but they emit 
high levels of carbon dioxide (CO,), nitrogen oxides (NO,) 
and sulfur oxides (SO,). Studies show that emissions of CO,, 
NO, and SO, from ships correspond to about 2, 11 and 4% of 
the global anthropogenic emissions respectively (Skjflsvik et 
al., 2000). Carbon dioxide is not toxic, but it contributes to the 
greenhouse effect (global warming). It is the main greenhouse 
gas emitted by ships. Nitrogen oxides contribute significantly 
to the ozone layer depletion2. They also form acids, which 
fall to earth as rain, fog, snow or dry particles, having harmful 
impacts on plants, animals and sea. Sulfur oxides, which are 
formed because of the high sulfur content of marine fuels (in 
order to reduce costs, ships generally use low quality fuels 
which tend to have a high sulfur content), are the major source 
of acid rain. To a lesser extent, other emissions from ships are 
carbon monoxide (CO), which is toxic to humans and animals; 
hydrocarbons, which contribute to the greenhouse effect and 
particulates, which can induce cardiopulmonary diseases and 
lung cancer. In harbour cities, emissions from ships constitute 
an important source of pollution. Furthermore, some gases 


may be transported in the atmosphere over several hundred of 
kilometers, contributing to air quality problems on land. 

Merchant ships in international traffic are subjected to 
International Maritime Organization (IMO). In particular, the 
IMO Marpol Annex VI (Regulations for the Prevention of 
Air Pollution from Ships) limits SO, and NO, emissions from 
marine engines. As SO, emissions are a function of the sulfur 
content of fuel, the IMO has limited the sulfur content in fuels. 
Concerning the NO,, the IMO has established a curve which 
indicates the maximum allowable NO, emission levels related 
to the engine speed. Figure | shows the maximum NO, emission 
levels for marine diesel engines built after 2000 (IMO Tier I), 
2011 (IMO Tier II) and 2016 (IMO Tier III, only in designated 
special areas called emission control areas). Apart from these 
international limitations, some regions have developed regional 
and national tougher limits. 

Due to this legislation, the reduction of NO, emissions is 
very important for the new marine engine designs. Several 
studies about NO, emissions and other pollutants from marine 
engines can be found in the literature. For example, those of Li 
et al. [1], Lü et al. [2], Holub et al. [3], Jayaram et al. [4], Yan et 
al. [5], Uriondo et al. [6], Bernecic and Кадопја [7], Maiboom 
and Tauzia [8], Duran et al. [9] and Millo et al. [10]. These 


D The greenhouse effect is a process by which thermal radiation from the planetary surface is absorbed by atmospheric greenhouse gases, 
which re-radiate part of this heat towards the surface again. As a result, the earth's temperature is higher than it would be if direct 


heating by solar radiation were the only warming mechanism. 


2 The ozone layer is a deep layer in the stratosphere. It is composed by molecules of ozone, O,, which filter out a lot of the sun's UV-B 
rays. Certain chemical substances can destroy some of the ozone, and the consequences are potentially damage to the life forms on 


Earth. 
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Fig. 1. IMO maximum allowable NO, emission for marine engines 
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and other investigations have provided a great progress in the 
marine emission control technologies. However, the design of 
an efficient marine engine requires a complete understanding 
of the combustion process, which has a noticeably influence 
on the emissions. There are many possible combinations 
between the operating parameters (compression ratio, injection 
characteristics, combustion chamber geometry, etc), and the 
effects on the engine performance and emissions produced 
by variations of these parameters are non linear and often 
opposed. In this regard, CFD (Computational Fluid Dynamics) 
offers an important tool to study the fluid flow in detail. CFD 
is a branch of fluid mechanics based on the subdivision of 
the computational domain into small elements overlaying the 
whole domain. For each subdivision, the governing equations 
are solved using numerical methods. 

This work proposes a CFD model to study the combustion 
process in the Wärtsilä 6L 46 four-stroke marine engine. The 
temperatures, pressures, advancing of the flame front and 
composition of the exhaust gas were obtained. The numerical 
model was validated with experimental measurements 
performed on a Wártsilà 6L 46 installed on a tuna fishing vessel, 
obtaining a good agreement. This work is a continuation of 
Lamas et al. [11], in which the exhaust, intake and compression 
strokes of the cylinder operation were simulated. 


CASE STUDIED 
Technical specifications 


The Wärtsilä 6L 46 studied in the present work is shown 
in Fig. 2. This is a four-stroke marine diesel engine with 
6 cylinders in line, 46 cm bore, 58 cm stroke and 96400 cm3 
cylinder displacement volume. Each cylinder has two intake 
valves and two exhaust valves, and the fuel injector is placed 
at the center of the cylinder head. Figure 3 outlines a cross- 
section of the engine. 

The present work was developed under the conditions 
studied in Lamas ef al. [11]. These are 96% load, 499.6 rpm 
and employing heavy fuel oil (RMG 380 according to ISO 
8217). Under these conditions, the power was 5212 kW and the 
consumption 172 g/kWh (9.96 g of fuel injected per cycle and 
cylinder). The in-cylinder gauge pressure is shown in Fig. 4, 
obtained using the engine performance analyzer MALIN 6000 
(Malin Instruments, Ltd.). 


Performance and combustion process 


The performance of this engine was described in Lamas et 
al. [11], where it was explained that combustion takes place 
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Fig. 2. Engine studied in the present work 
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Fig. 3. Cross section. Adapted from Wärtsilä [12] 
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Fig. 4. In-cylinder pressure experimentally measured 


towards the end of the compression stroke, with the piston 
near top dead center (TDC). As this is a direct injection engine, 
the fuel 1s injected directly into the combustion chamber. The 
combustion chamber is a shallow bowl in the piston crown 
shown in Fig. 5, which illustrates a photograph ofthree pistons 
of this engine. 


Fig. 5. Pistons of the engine studied in the present work 


The injection process is shown schematically in Fig. 6(a). 
The fuel is introduced by a 10 holes injector, Fig. 6(b), in the 
form of a spray of liquid droplets. As the air contained in the 
cylinder is at high pressure and temperature, the fuel vaporizes 
and mixes with this air. After that, ignition takes place. 
Atomization, vaporization, fuel-air mixing and combustion 
continue until all the necessary fuel has passed through each 
process. In addition, mixing ofthe air remaining in the cylinder 
with burning and already burned gas continues throughout the 
combustion and expansion processes. 


Injector 


Piston 


Cylinder 


b) 


Fig. 6. a) Schematic representation of the spray injection; 
b) photograph of the injector 


COMBUSTION CHEMISTRY AND 
GOVERNING EQUATIONS 


Combustion is a very complex phenomenon which involves 
different disciplines such as fluid mechanics, chemistry, 
thermodynamics, heat transfer, etc. It is defined as an oxidation 
process of a fuel (in this case diesel, C,,H,,) and an oxidant 
(in this case oxygen, O,), converting chemical energy into 
heat energy. 


Ifa fuel is burnt stoichiometrically and ideally, this would 
be oxidized with air (which has a molecule of oxygen per 3.76 
molecules of nitrogen) by means of the following reaction: 


4C, H,, +7100, +3.76N,) > 
— 48CO, +46Н,0 +71-3.76N, 


This reaction 1s very simple because it only involves 
five chemical species, C,,H,,, О,, №, CO, and H,O and four 
elements, C, H, O and N. Unfortunately, this ideal mechanism 
never occurs and real combustion processes are governed by 
hundreds of reactions and tens of chemical species. Currently, 
there are no computational resources to model all the reactions 
and species involved in the combustion process. The solution 
is to simplify the physical phenomenon which governs the 
chemical reactions, explained in what follows. 


(1) 


Governing equations 


In this engine, the cylinder is full of air before the injection 
of fuel, and combustion takes place as the fuel is being injected. 
This is known as non-premixed combustion (in contrast to 
premixed combustion, in which fuel and oxidant are mixed 
before combustion takes place). In non-premixed combustion, 
the rate of combustion is controlled by the rate at which fuel 
and air mix. Under such hypothesis, the equation which 
characterizes the propagation of the flame front is given by: 


д 0 | u, ôf 
==. |: 
д Pa) А 


1 


д 
Lipi 2 
XD (2) 


а, дх, 

where p is the density, o; the turbulent Schmidt number and f 
the mixture fraction. This is the local mass fraction of burnt and 
unburnt fuel stream elements (C, H, etc) in all the species (CO,, 
H,O, О,, etc). For the following chemical reaction, in which C 
represents the fuel, O the oxidant and P the products: 


cC +00 — (с + o)P (3) 
the mixture fraction is given by: 
f=c/(c + o) (4) 


In addition to this equation, the equations of conservation 
of mass, momentum and energy are also necessary, Eqs. (5) 
to (7) respectively: 
др @ 
ob (pu) =0 (5) 
Ot Ox 


0 0 
^ (рц) + ax, (pu,u,)= 
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In the equations above, p is the pressure, t the viscous 
stress tensor, H the enthalpy, и the, o the Prandtl number, u, 
the turbulent viscosity, o, the turbulent Prandtl number and $, 
a source term to include the radiation heat transfer. The term 
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—pu;u, represents the Reynolds stresses, which were computed 
by means of the К-е turbulence model. 

Concerning the characterization of the chemical species 
involved, these were computed using the concept of chemical 
equilibrium (the chemical kinetics is so fast that equilibrium 
18 reached quickly). The following 20 chemical species were 
considered: C,;H,, (fuel), О,, Nj, СО, Н.О, CO, С, CH,, O, H, 
H,, Н,О,, N, OH, HO,, HNO, HONO, C,H, HCO and CHO. 

NO, is formed due to the high temperatures reached during 
the combustion process. Around 1500°C and above, enough 
thermal energy is available to dissociate oxygen and nitrogen, 
which recombine to form NO,, Versteeg and Malalasekera [15]. 
NO, emissions from ships are relatively high because most 
marine engines operate at high temperatures and pressures 
without effective NO, reduction technologies. The reactions 
of formation of NO, are not fast enough to consider chemical 
equilibrium. For this reason, its treatment was decoupled from 
the combustion model. Among the NO, components, NO is the 
main pollutant (nearly 100% of NO, is NO). NO is generated 
by the reactions proposed by Zeldovich [13], reactions (8) and 
(9), and Lavoie [14], reaction (10). 


К, 
O+N, > NO+N (8) 


N+0, & NO+0 (9) 


N+OH & NO+H 10 


The rate of formation of NO is given by the following 
expression: 


d[NO]/dt = [ОЈ] + ҸО] 
+ k[N][OH] - к ,[NO][N] + 


~ КОО] - k ,[NO][H] 


where [NO], [О], [Ч], [N], [0,], [OH] and [H] are the 
concentrations, mol/m:. k,, k,, K, are the rate constants and k, 
k, Kk, are the corresponding reverse rates. The forward and 
backward reaction rate constants for these three key reactions 
are given by: 


(11) 


k - 1 8 xl oua ew 
1 . 
k, =1.8x10 Тес 699m 


бе 7.110760 


k =3.8х10%еС 95D (12) 
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As can be seen, the reactions above are very dependent on 
the temperature. The NO formation rate is much slower than 
the combustion rate and most of the NO is formed after the 
completion of the combustion, due to the high temperatures 
present in the combustion zone. 


NUMERICAL PROCEDURE 


The numerical computations were performed by means of 
the commercial software ANSYS Fluent. It was carried out 
a period of time from 30? before TDC to the exhaust valves 
opening, 127? after ТОС. АП valves remained closed during 
the entire simulation. The results of Lamas et al. [11] were 
continued to 30? before TDC and employedas initial conditions 
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of the present work. Concerning the boundary conditions, the 
heat transfer from the cylinder to the cooling water was modeled 
as a convection type: 


q= h(T pas m Tater) 


where q is the heat transferred, T,,, is the in-cylinder 
temperature, T,,,,., is the cooling water temperature (78°C) 
and h is the heat transfer coefficient, given by the following 


expression, Taylor [16]: 
h = 10.4 kb-!^ (u (14) 


where b is the cylinder bore, k the thermal conductivity of the 
gas, у... the mean piston speed and v the kinematic viscosity 
of the gas. Substituting values into the above equation yields 
h=4151 W/mK. 


13) 


йыда” yp 


In order to model the movement of the piston, a moving 
grid was employed. Figure 7 shows the grid at the start of the 
simulation, 30? before TDC, and at TDC. As can be seen, only 
the combustion chamber was simulated. 


-30° TDC 


Fig. 7. Computational grid 
RESULTS AND DISCUSSION 


The spray and mass fraction of fuel are shown in Fig. 8. 
Under the conditions studied, the injection of fuel takes place 
from -12° to -1° crankshaft angles. Figure 8 represents - 10°, -5° 
and 0° (PMS). As can be seen, the fuel spray enters the cylinder 
as liquid drops (which were simulated as 5 um diameter). These 
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Fig. 8. Spray and mass fraction of fuel 


drops are represented as black points in the figure. The structure 
of each fuel spray is that of a narrow liquid-containing core 
surrounded by a much larger gaseous-jet region containing 
fuel vapour. The flame spreads rapidly as the fuel is injected, 
producing a very fuel rich zone close to the fuel пр (-10° 
and -5?). After the injection of fuel has finished, the flame is 
extinguished progressively (0?). 

The temperature field is shown in Fig. 9. The fuel is injected 
at 127?C, creating a cold temperature in the zone closed to the 
injector. In the rest ofthe cylinder, the temperature increases as 
the crankshaft angle is advanced due to the combustion progress 
and the compression supplied by the piston. 


0° (PMS) 
[ ^ 
2 Y 
Fig. 9. Temperature field (°С) 


Figure 10 illustrates the comparison of simulated and 
measured in-cylinder pressures. A satisfactory agreement 
between the experimental and numerical results was 
obtained. 

200 


150 
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Pressure [bar] 


-30 0 30 60 90 120 
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Fig. 10. In-cylinder pressure numerically and experimentally obtained 


The emissions of carbon dioxide, carbon monoxide, 
nitrogen oxides and hydrocarbons are shown in Table 1. The 
NO, emissions were measured experimentally using the G4100 
analyzer (Green Instruments). According to this apparatus, 


the engine emits 1115 ppm of NO,, which suppose an error of 
13.8% with respect to the numerical result given in Table 1. This 
error is quite acceptable considering the difficulty of simulating 
the combustion process. As the CFD code provides the quantity 
of mass inside the cylinder, the NO, emissions in g/kWh can be 
easily calculated, obtaining a value of 11.8. This value sticks 
the IMO TIER I but not IMO TIER II, applicable to engines 
mounted in vessels built after 1 January 2011. Several solutions 
to decrease the NO, levels in this engine were proposed by 
Millo et al. [10]. 


Tab. 1. Emissions from the engine numerically obtained 


CO, CO 
(96 vol.) | (ppm* vol.) 


NO HC 


(ppm vol.) (ppm vol.) 


961 569 


*ppm: parts per million 


It is interesting to study the NO, emissions against the 
start of injection (SOI). For the purpose, the SOI was varied 
between -13? and -10?. The experimental and numerical 
results are shown in Fig. 11. The peak pressures numerically 
and experimentally obtained are also shown in the figure. As 
can be seen, delaying the start of injection reduces the NO, 
emissions and the peak pressure, both for the experimental and 
numerical results. The NO, emissions are very sensitive to the 
temperature. The decrease in NOx emissions with injection 
retard is related to the decrease in combustion pressure and 
thus temperature. Unfortunately, the power 1s lower as the 
pressure decreases. 
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Fig. 11. NOx emission varying with injection timing 


CONCLUSIONS AND FURTHER 
DEVELOPMENTS 


-11.0 -10.5 


The present paper proposes a CFD model to simulate the 
combustion process in the Wärtsilä 6L 46 marine engine. The 
strongest motivation is given by the current legislation, for 
which the most important gas component that must be reduced 
in exhaust emissions from marine engines is NO,. 

The pressure and temperature fields were obtained, as well 
as the composition of the exhaust gas. In order to validate 
this work, numerical results were satisfactory compared with 
experimental ones performed on a Wärtsilä 6L 46 installed 
on a tuna fishing vessel. This work provides a model to study 
several parameters such as injection pressure, compression 
ratio, injection timing, etc. In future works, the purpose is to 
employ this model to design low emission engines, analyzing 
measures such as water addition, retarded injection, modulated 
injection, exhaust gas recirculation, etc. 
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NOMENCLATURE 


— Cylinder bore 

=  Mixture fraction 

— A Heat transfer coefficient 
Enthalpy 

— Thermal conductivity 

— Pressure 

— A Heat transfer 

— Time 

= Velocity 


Е бою ә ъс 
| 


Special characters 


u = Turbulent viscosity 

v = Kinematic viscosity 

p - Density 

o, —  Turbulent Prandtl number 

с; = Turbulent Schmidt number 

Tj — Stress tensor 

Subscripts 

i = Cartesian coordinate (i= 1, 2, 3) 
j - Cartesian coordinate (j = 1, 2, 3) 
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ABSTRACT 


In this paper are characterized problems of high strength steel weldability in underwater wet welding 

conditions. Water as a welding environment intensifies action of unfavourable factors which influence 

susceptibility to cold cracking of welded steel joints. The susceptibility to cold cracking of $355J2G3 steel 

and S500M steel in wet conditions was experimentally estimated (by using Tekken test). It was concluded 

that the steels in question are characterized by a high susceptibility to formation of cracks in welds. 

Usefulness of the proposed Temper Bead Welding technique (TBW) was experimentally verified as a method 
for improving weldability of the steels in the analyzed conditions. 


Key words: underwater welding; steel weldability, cold cracking; wet welding; cracks in welds; improving weldability 


INTRODUCTION 


According to the most general classification, the underwater 
welding processes can be divided into dry and wet ones [1-5]. 
The dry methods consist in local isolation of welding area 
from surrounding water by means of special chambers in 
which atmospheric pressure (isobaric welding) or an elevated 
pressure resulting from welding depth (hyperbaric welding) 
is present. A characteristic feature of the wet welding is 
direct underwater contact of diver-welder, electrode and base 
material. An intermediate variant between wet and dry welding 
is application of a local dry chamber (dry spot welding) [1-3]. 
Applicability of particular processes in underwater welding 
conditions 15 characterized in Tab. 1. 


Among wet welding methods, the manual metal arc welding 
is most often used, hence majority of steel weldability test 
results accessible in the subject-matter literature concern the 
above mentioned process [4-6, 11-20]. 

Transferring the welding process to water environment 
results first of all in increasing diffusible hydrogen content in 
deposited metal, as well as in increased cooling rate. Manual 
metal arc welding in wet conditions generates the diffusible 
hydrogen content in deposited metal, ofthe order ofa few dozen 
ml/100g Fe, depending on a type of shielding and welding 
conditions [14, 15]. The tests [15, 16] performed with the use 
of rutile-shielded electrodes, the most often applied in such 
conditions, showed that hydrogen content in deposited metal 
did not depend on a wetting degree of the shielding, but heat 
input was the decisive factor (Fig. 1). 


Tab. 1. Possibility of application of welding processes in water environment [2-10] 


No. Process 


- 


ММА (111) уеѕ yes 


N 


Dry welding 


І Local chamber 
Wet welding method welding 


possible 


SAW (121) yes yes 


MIG/MAG (131/135) yes 


possible yes 


FCAW (136/138) yes 


yes yes 


FCAW-S (114) 


possible yes 


possible 


Plasma arc welding (15) 


possible no 


possible 


Laser welding (521) 


possible no yes 


Friction welding (42) yes 


yes no 


3 
4 
5 
6 
7 
8 
9 


Explosion welding (441) no 


yes no 


= 
i=) 


Stud arc welding (78) 


possible yes 


possible 
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Fig. 1. Effect of heat input of ММА wet welding on diffusible hydrogen 
content in deposited metal [16] 


Cooling rates of joints made underwater are much shorter 
compared with those for in - air welding conditions. The 
increase of cooling rate contributes to forming brittle structures 
in heat affected zone (HAZ) and to increasing values of residual 
stresses. Such conditions result in that welded structural steel 
joints made directly in water are characterized by a high 
susceptibility to forming cold cracks. The cracks appear even 
in joints made of steel of a low carbon equivalent Ce < 0.3 96. 
In the HAZ of such steel joints the increase of hardness up to 
the level of 350 + 400 HV, as well as many hydrogen micro- 
cracks were observed [13, 14]. In Fig. 2 are presented the results 
of Tekken tests performed underwater, which were used for 
development of an equation for determining the minimum heat 
input value which protects against cold cracks forming [4]. 


20 


16 


Heat input q [kJ/cm] 


12 


0.35 040 045 050 0.55 060 065 070 0.75 
Cracking parameter (P-T/600) [% 


Fig. 2. Results of Tekken tests, underwater wet MMA welding: 
о — no cracks, х — internal cracks, € — internal and external cracks [4] 


The equation is of the following form: 


q = 200 х (P — T/600) — 67 


where: 
а - heat input, [kJ/cm], 
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P = Се + H,/60 + R/200000, 


T - initial temperature, [9C], 

Ce — carbon equivalent acc. International Institute of Welding, 
[76], 

Hp = diffusible hydrogen content in deposited metal, [ml/100g 
Fe], 


R - restraint intensity, [MPa]. 


The tests [17] on Implant specimens, performed in the 
conditions similar to those for Tekken test, showed that critical 
stresses were higher than yield point ofthe tested steel. It goes to 
show that the steel was crack-resistant, which 15 not consistent 
with the results given in [4]. 

A great number of variable factors which act during 
underwater welding makes it difficult to formulate safe 
welding conditions, i.e. those possible to protect joints against 
cracking. However on the basis of many tests the following 
approximate criteria for safe wet welding of non-alloy steels, 
were determined [18]: 

- base material carbon equivalent: Ce < 0.30 95, 
- heat input: greater than 1.5 kJ/mm. 


In the case of application of austenitic electrodes value of 
welded steel carbon equivalent may be even greater — up to 
0.40 95. In this case hot cracks should be expected to occur 
in welds [4, 19]. Increase of toughness of weld metal at wet 
welding conditions can be also achieved by applying low-alloy 
electrodes containing nickel. The optimal Ni content in weld 
is considered as 2% [19]. 

Another limitation method of susceptibility to forming 
cold cracks is application of temper bead welding technique 
(TBW) [20, 21]. In the TBW a reduction of post-welding 
residual stresses as well as a decrease of hardness in HAZ is 
achieved. The technique is applied to thick sections, i.e. where 
it 18 not possible to perform heat treatment operation after 
welding. It consists in overlaying successive beads onto weld, 
that introduces local heat treatment in the area of prior made 
weld layer. In HAZ ofthe layer structural transformations may 
occur in the temperature range higher than A, temperature. 
ASME rules require to make six layers, and for three first layers 
increased value of heat input is often used [21]. In the case of 
implementation of TBW technique for preventing against cold 
cracking it is important to determine a time interval between 
instants of overlaying a tempered bead and tempering one. 
Depending on a type of applied electrode the time up to crack 
forming at wet welding may reach from 3 min to 2 hours. 
Successive parameters important from the point of view of 
effectiveness of the process are: value of heat input in every 
bead as well as distance between bead axes (pitch) [21]. 


EXPERIMENTAL 


The tests were aimed at determination of susceptibility to 
forming cold cracks in joints of high strength steels, made by 
using covered electrodes in wet welding process. The tests were 
performed in accordance with the following scheme: 

- preparation of welding stands in air and underwater 
environment 

- preparation of Tekken tests 

- welding of test joints 

- visual testing and penetrant testing 

- macroscopic and microscopic metallographic 
examinations 

- Vickers hardness measurements. 


Additionally, to assess usefulness of application of TBW 
technique for improving steel weldability in water environment, 


Tab. 2. Chemical composition of S355J2G3 and S500M steel, wt 96 


$355J2G3 


S500M 


test padding welds were made and subjected to macroscopic 
metallographic examinations and hardness measurements. 

For the tests 15 mm plates of S355J2G3 and S500M steel 
were used. Their chemical compositions and mechanical 
properties are given in Tab. 2 and 3. 


Tab. 3. Mechanical properties of $355J2G3 and S500M steel 
Yield point | Tensile strength | Elongation 
R, [MPa] R, [MPa] А, [v6] 
8355J2G3 464 5T] 22.1 
S500M 525 619 20.5 


The test welds were made by using Lincoln Electric OMNIA 
(E 42 ORC 11) electrodes of 4 mm diameter. They are general 
application rutile electrodes for welding in all positions, of the 
chemical composition: C — 0.07%, Mn — 0.5%, S1— 0.5%, and 
the mechanical properties: R, — min. 420 MPa, К, — 500 + 640 
MPa, А; — min 20%. 

The Tekken test joints of S500M steel were made in water 
and air environment at various heat input values. The Tekken 
test joints of S355J2G3 steel were made only at underwater 
conditions. The test joints were prepared in compliance with 
the standard guidelines [22], on the underwater welding stand 
at low water depth [3]. In Tab. 4 are presented the conditions 
for performing Tekken tests on selected joints. 


Penetrant testing 

The visual and penetrant testing [23-25] of the test 
joints were performed after 72 h from finishing the welding. 
Occurrence of cracks in welds of all the specimens made 
in water environment was observed, but in S355J2G3 steel 
joints the cracks were placed close to fusion line, whereas 
in the specimens made of S500M steel — in weld axis. In the 
Tekken specimens welded in air environment no cracks were 
revealed. 


Macroscopic metallographic examinations 

The macroscopic metallographic examinations were 
performed in accordance with the standard guidelines [26]. 
Example photos of the cross-sections are presented in Fig. 
3-5. In the specimens made underwater extensive cracks and 


a porosity were sometimes found. Tekken specimens of S500M 
steel welded in air environment were characterized by a much 
lower susceptibility to cracking, however in one of them some 
cracks were observed. The results obtained from the Tekken 
tests are presented in Tab. 5. 


Fig. 3. Cross-section of the W1 (S500M) test joint with the visible crack 
initiated in weld root 


Fig. 4. Cross section of the 1A (S355J2G3) test joint with the visible crack 
initiated in fusion line 


Fig. 5. Cross-section of the P12 (S500M) test joint with the visible crack 
initiated in weld root 


Tab. 4. Conditions of Tekken tests 


Environ 


No. of 
ment 


parameters 


Welding Welding 


(түз Heat input value 


specimen 


5500М 


UIV] 


q [kJ/mm] 


I [A] t [s] 


S500M 


S500M 


S500M 


S500M 


535517263 


53551203 
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No. of 
specimen 


Environment 


Penetrant 
testing 


crack 


Tab. 5. Results of Tekken tests 


Macroscopic 
metallographic 
examinations 


crack 


Hardness 
measurements 


Percentage of cracks 


HV, max 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


no crack 


no crack 


no crack 


no crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


crack 


Fig. 6. Microstructure of S500M and 853557283 steel joints welded underwater; a) base material, b) weld, c) overheated zone in HAZ. Magnification x 200 


70 
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crack 


Microscopic metallographic examinations 

In Fig. 6 are presented typical structures of base material, 
weld and HAZ at fusion line for the tested joints made 
underwater. The S500M steel structure consists of ferrite grains 
of very different sizes with a trace of pearlite content whereas 
53557203 steel shows a ferritic-pearlitic structure of distinct 
banding. The weld structure of both the steels 15 similar and 
consists of ferrite of column arrangement with an outline of 
Widmanstátten structure. In the HAZ of S500M steel, quasi- 
pearlite structures with acicular ferrite bands reaching down 
the grains, are observed, whereas in the HAZ of 53551203 
steel acicular quench structures are visible. 


Hardness measurements of welded joints 

The hardness measurements were performed according to 
the standard [27] under 98 N load (HV10). Example results are 
shown in Fig. 7 through 9. In the S500M steel joints made both 
in the air and underwater, hardness values for weld and base 
material were close to each other and equal to about 200 HV10. 
In HAZ only a slight increase in hardness up to 240 HV10 was 
observed for air environment welding conditions, whereas 
the transferring of the welding process to water environment 
resulted in further increase of maximum hardness of HAZ up 
to about 280 HV10. 

However the maximum hardness values in HAZ of the 
8355J2G3 steel joints were significantly greater. A higher 
carbon content in this steel in association with an increased 
cooling rate of the joints resulted in the quenching of HAZ 
material and increase of its hardness up to over 400 HV10. 

300 
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Fig. 7. Hardness distribution across P6.2 test joint 
of S500M steel welded in the air, HVmax — 242 
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Fig. 8. Hardness distribution across W2.1 test joint 
of 5500М steel welded underwater, HVmax = 283 
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Fig. 9. Hardness distribution across 1A.1 test joint 
of S355J2G3 steel welded underwater, HVmax = 425 


Application of temper bead welding (TBW) 

A preliminary assessment of usefulness of application of 
the tempering bead technique to improve weldability in water 
environment was performed by using $355J2G3 steel plates. 
On the testing plate of 12 x 100 x 200 mm dimensions padding 
welds were made with the use of Omnia electrodes at 10 min 
interval and 4096 pitch. Cross-sections of the padding welds 
were subjected to the macroscopic metallographic examinations 
and hardness measurements. Parameters of the padding welding 
are given in Tab. 6. A photo of polished macroscopic cross- 
section of test beads with indication of padding weld laying 
sequence is shown in Fig. 10. 


Tab. 6. Welding parameters for $355J2G3 steel 


Welding 
current 
I [A] 


Arc 
voltage 
U [V] 


Welding 
speed v 
[mm/min] 


Heat input q 
[kJ/mm] 


Fig. 10. Cross section of the test weld beads; S355J2G3steel; 40% pitch 


In Fig. 11. is shown location and designation of areas in 
which hardness measurements were made, and in Tab. 7 and Fig. 
12 are presented results of hardness measurements performed 
at 2 mm distance measured from specimen surface. 

The presented results of hardness measurements for padding 
welds make it possible to conclude that the temper bead welding 
in underwater welding allows to lower maximum hardness 
in HAZ to a safe level as regards cold crack forming. In Fig. 
12 two characteristic points of increased hardness (swcap 
and swcbp) can be observed. They appear at fusion line and 
simultaneously very close to edges of the padded elements. 
In this place intensive heat sinking occurs, a.o. due to bubble 
boiling, that leads to forming high hardness structures. In 
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Tab. 7. Results of hardness measurements of S355J2G3 steel test plates 


|. | mr  swean| sweap | mna | wsweb | mnb | swebp | swebn 


[vests [sse [тулеш [rap 


ную | эз | xs | 23 | 


mna wswcb mnb 


Fig. 11. Location and designation of tested areas on testing plate: 
mr — base material, mna — 1* padding weld, mnb — 2"! padding weld, 
swca – HAZ ој l“ padding weld, swcb —HAZ of 2" padding weld, 
wswca — area of overlapping swca and swcb, wswcb — area of swcb 
overlapping mna, n — normalization area, p — overheated area, 
А-А – hardness measurement line 
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Fig. 12. Hardness distribution along А-А line (see Fig. 11) on S355J2G3 
steel test plate. The HVmax ој the tempered zone in HAZ ој I* weld bead 
(swca2p) is marked by red horizontal line. 
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making multilayer joints such phenomenon can be minimized 
by providing heat again during laying successive tempering 
beads — in the case in question — it means making the last bead 
between two padding welds. 


DISCUSSION 


In all the test joints made underwater, cracks in welds were 
observed. The phenomenon is in compliance with information 
contained in the literature sources concerning weldability of the 
tested steels [28]. Specimens made in the air were characterized 
by a much lower susceptibility to crack forming. Crack 
occurrence in the joints made in the air resulted from using rutile 
electrodes. Application of such electrodes was conditioned 
in association with real underwater welding conditions with 
the use of the MMA process as well as convergence with 
requirements of the PN-EN ISO 15618-1 standard concerning 
certification of divers-welders for wet welding [29]. In the 
above mention standard, one of the recommended welding 
electrodes 15 that with rutile coating. Welding with rutile 
electrodes both in the air and underwater is a process of a high 
hydrogen content. Presence of hydrogen in steel weldments 
leads to degradation oftheir structure and mechanical properties 
[30]. On the basis ofthe previous tests the diffusible hydrogen 
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content in deposited material obtained in air environment can 
be estimated to be about 35 ml/100g Fe, whereas in the case 
of underwater welding — 50 ml/100 g Fe [15, 31]. 

In the wet welding conditions for S500M steel the increase 
of HAZ hardness from 240 + 250 HV10 (characteristic for 
welding in the air) up to 260 + 300 HV10, took place. Maximum 
hardness values in HAZ of S500M steel are lower than 350 
HV10 value is considered as critical from the point of view 
of susceptibility to cold cracking. However the test joints 
of S355J2G3 steel were characterized by the much greater 
maximum hardness in HAZ of values over 400 HV10. It results 
from their welding conditions (a. o. greater cooling rates - 
consequently - shorter cooling times t,,;), which, at the carbon 
equivalent Ce = 0.44% for S355J2G3 steel, led to forming 
quenched structures. For Tekken specimens of 8355J2G3 
steel the initiating of cracks in fusion line 18 characteristic 
(Fig. 4), which is confirmed by presence of brittle structures in 
this area. Hardness distribution over cross-section of Tekken 
joints made of S500M steel complies with acceptance criteria 
used for approving welding technologies in the air. Limitation 
of maximum hardness in HAZ down below 350 HV10 value 
results from low carbon content in the steel as well as carbon 
equivalent Ce — 0.3096. The presented results indicate that 
low-carbon steels after thermo-mechanical processing is 
applicable for marine and off-shore engineering structures 
when underwater welding is required. 

The application of the TBW technique to $355J2G3 steel 
allowed to reach satisfactory results, that is manifested by 
the hardness measurement results in trouble areas of padding 
welds (Tab. 7). The technique in question is greatly useful 
in welding elements made of a steel which requires a heat 
treatment after welding. In wet welding conditions or dry 
hyperbaric ones such treatment would be very expensive and 
in most cases impossible to perform. TBW technique makes it 
possible to obtain, by applying relatively simple technological 
operations, welded joints which satisfy criteria of subject- 
matter standards or requirements of classification societies. It 
may find application e.g. to repair work of underwater piping 
lines. It should be added that effectiveness of application 
of tempering beads depends on many factors e.g. welding 
position, welding parameters (heat input), arc length, torch 
setting angle, interval between padding weld laying and weld 
pitch. However further investigations are required to develop 
TBW technology for wet welded joints. The successive tests 
should take into account factors, that would make it possible to 
determine optimum conditions for application of temper bead 
technique to underwater welding process. 


CONCLUSIONS 


1. The joints of S355J2G3 and S500M steel, made underwater 
in the restraint conditions showed high susceptibility to 
forming cold cracks in welds. 

2. Inall the Tekken test joints of S500M steel, made with the 
use of wet welding method cold cracks were revealed, and 
in the joints made in the air such cracks were found in only 
one of the two test joints. 


3. Maximum hardness in HAZ of the S500M steel joints did 
not exceed 300 HV, whereas that in HAZ of the S355J2G3 
steel joints was over 400 HV10 in value. 

4. It was proved that the TBW technique was effective for 
the joints of high strength steel welded underwater. The 
application of the TBW resulted in the decreasing of 
maximum hardness in HAZ of 83551203 steel joints down 
to the value below 350 HV10. 


BIBLIOGRAPHY 


1. The Standard AWS D3.6M:2010 Underwater Welding Code. 

2. Labanowski J.: Development of underwater welding techniques. 
Przeglad Spawalnictwa (Welding Technology Review), 10/2008 
(in Polish). 

3. Labanowski J., Fydrych D., Rogalski G.: Underwater Welding 
— a review. Advances in Materials Science, 3/2008. 

4. Christensen N.: The metallurgy of underwater welding. 
Proceedings of the International Conference ,,Underwater 
Welding", Trondheim, Norway 1983. 

5. Cotton H. C.: Welding under water and in the splash zone 
— a review. Proceedings of the International Conference 
„Underwater Welding", Trondheim, Norway 1983. 

6. Ibarra S., Grubbs C. E., Liu S.: State of the art and practice 
of underwater wet welding of steel. International workshop on 
underwater welding of marine structures. New Orleans, USA, 
1994. 

7. Zhang X., Ashida E., Shono S., Matsuda F.: Effect of shielding 
conditions of local dry cavity on weld quality in underwater 
Nd:YAG laser welding. Journal of Materials Processing and 
Technology. Vol. 174, Issues 1-3, May 2006. 

8. Fydrych D., Rogalski G.: Effect of underwater local cavity 
welding method conditions on diffusible hydrogen content in 
deposited metal. Welding International Vol. 27, Issue 3, March 
2013, pages 196-202. 

9. Ambroziak A., Gul B.: Underwater friction bonding of overlap 
joints by means of a plasticized steel pin (FHPP). Przeglad 
Spawalnictwa (Welding Technology Review), 9-10/2006 (in 
Polish). 

10.Ambroziak A., Gul B.: Investigations of underwater FHPP for 
welding steel overlap joints. Archives of Civil and Mechanical 
Engineering 2/2007. 

11.Rodriguez-Sanchez J. E., Rodriguez-Castellanos A., Perez- 
Guerrero F., Carbajal-Romero M. F., Liu S.: Offshore fatigue 
crack repair by grinding and wet welding. Fatigue and Fracture 
of Engineering Materials and Structures, 34/2010. 

12.Pessoa E., Bracarense A., Zica E., Liu S., Perez-Guerrero F.: 
Porosity variation along multi-pass underwater wet welds and 
its influence on mechanical properties. Journal of Materials 
Processing Technology. Vol. 179, Issues 1—3, 20 October 2006, 
pp. 239-243. 

13.Bohme D., Eisenbeis C.: Jnvestigation into the credibility of the 
implant test when used to assess the cold cracking sensitivity 
of underwater wet welds. Proceedings of the International 
Conference ,,Welding Under Extreme Conditions”, Helsinki, 
Finland 1989. 

14.Hoffmeister H., Kuster K.: Process variables and properties 
of underwater wet shielded metal arc laboratory welds. 
Proceedings of the International Conference „Underwater 
Welding", Trondheim, Norway 1983. 

15.Fydrych D., Rogalski G.: Effect of shielded-electrode wet 
welding conditions on diffussion hydrogen content in deposited 
metal. Welding International Vol. 25, Iss. 3, March 2011, p. 166- 
171. 


16.Suga Y.: Effect of diffusible hydrogen on mechanical properties 
of underwater welded joints — study on improving the 
mechanical properties of underwater welded joints (the 17 
report). Transactions of the Japan Welding Society 10/1985. 

17.Brink S. H., Boltje С. W.: Cold cracking susceptibility of 
welds obtained by wet underwater welding. Proceedings of the 
International Conference „Underwater Welding", Trondheim, 
Norway 1983. 

18.Fydrych D.: Cold cracking of steel welded in water environment. 
Przeglad Spawalnictwa (Welding Technology Review), 10/2012 
(in Polish). 

19.Liu S.: Fundamentos de saldatura humeda. Seminario de 
saldatura humeda. Campeche, Mexico 1999. 

20.Global Divers & Contractors : Joint industry underwater 
welding development program. Phase I - final report. Colorado 
School of Mines, 1995. 

21.Lomozik M.: Morphology and toughness of heat affected zone 
regions of steel welded joints in the aspect of temper beads 
application. Uczelniane Wydawnictwa Naukowo-Dydaktyczne 
AGH (Scientific and Didactic Publishing House, Cracow Mining 
and Metallurgy Academy), Kraków 2007 (in Polish). 

22.The Standard PN-EN ISO 17642-2:2005. Destructive tests on 
welds in metallic materials. Cold cracking tests for weldments. 
Arc welding processes. Part 2: Self-restraint tests. (in Polish). 

23.The Standard PN-EN ISO 17637:2011. Non-destructive testing 
of welds. Visual testing of fusion-welded joints (in Polish). 

24.The Standard PN-EN 571-1:1999. Non-destructive testing. 
Penetrant testing. Part 1: General principles (in Polish). 

25.The Standard PN-EN ISO 23277:2010. Non-destructive testing 
of welds. Penetrant testing of welds. Acceptance levels (in 
Polish). 

26.The Standard PN-EN 1321:2000. Destructive test on welds in 
metallic materials. Macroscopic and microscopic examination 
of welds (in Polish). 

27.The Standard PN EN ISO 9015-1:2011 Destructive tests on 
welds in metallic materials. Hardness testing. Hardness test on 
arc welded joints (in Polish). 

28.Tasak E.: Weldability of steel (in Polish ). Wydawnictwo Fotobit 
(Fotobit Publishers), Cracow 2002. 

29.Rogalski G., Labanowski J.: Certification of divers-welders 
for underwater wet welding in hyperbaric conditions. Biuletyn 
Instytutu Spawalnictwa (Bulletin of Polish Welding Institute), 
1/2011 (in Polish). 

30.Cwiek J.: Hydrogen assisted cracking of high-strength weldable 
steels in sea-water. Journal of Materials Processing Technology, 
Vol. 164-165 (2005), p. 1007-1013. 

31.Fydrych D., Labanowski J.: Determining diffusible hydrogen 
amounts using the mercury method. Welding International Vol. 
26, Iss. 9, September 2012, p. 697-702. 


CONTACT WITH THE AUTHORS 


Dariusz Fydrych, Ph.D., 
Jerzy Labanowski, Assoc. Prof., 
Grzegorz Rogalski, Ph.D., 
Faculty of Mechanical Engineering 
Gdansk University of Technology 
Narutowicza 11/12 
80-233 Gdansk, POLAND 
e-mail: jlabanow@pg.gda.pl 


POLISH MARITIME RESEARCH, No 2/2013 73 


POLISH MARITIME RESEARCH 2(78) 2013 Vol 20; pp. 74-81 
10.2478/pomr-2013-0019 


An analysis of the welding-induced deformation 


of ship panels using a statistical tool 


Huggo S. Batista, M.Sc., 
Tadeusz Graczyk, Assoc. Prof., 
West Pomeranian University of Technology, Szczecin 


ABSTRACT 


This paper proposes a practical way to predict the welding-induced deformation of 5-mm-thick steel ship 
panels prefabricated in a Polish shipyard, developed by means of a number of experiments and a body of data 
collected on the butt welding line using the submerged arc welding technique. The program of the experiments 
was so designed as to test different welding energy conditions as well as pressure values generated by 
the devices used during the welding. Subsequently, different welding trials were carefully carried out 
and data collected to find the most important parameters which contribute to welding deformation. The 
data obtained were analysed using the Design-Expert® statistical software, and consideration was given 
to the data arrangement method that would allow for the most appropriate description of the welding 
deformation. As a result of the foregoing efforts, an approach to predicting the welding deformation, 
based on a numerical equation including the main welding parameters, was developed. The accuracy of 
the equation was partially evaluated and the results turned out to be satisfactory, as the actual and the 
predicted values were comparable. Another important advantage achieved was the arrangement of data 
used to compute the welding deformation. 


Keywords: ship production technology; welding deformation; design of experiment; statistics 


INTRODUCTION 


Analysing the contemporary ship production requires some 
attention on issues related to design of ships, the capacity of 
shipyards and work stations, as well as production technology 
due to the fact that business factor is related to them. In 
consequence these subjects are in permanent development 
because of their significance for the production costs besides, 
schedule and quality. Consequently, there is a need in the 
shipbuilding industry for improvements to shorten delays and 
avoiding rework due to quality problem or even creating new 
solutions for the production processes. 

The following issues are examples of possible improvements 
in the ship production technology: 

- development of new raw materials that can be used in the 
shipbuilding industry; 

- reduction of rework related to deformations occurring after 
welding (straightening processes); 

- reduction of the time required to prepare cut outs in steel 
parts; 

- development of technologies for fitting parts to be welded 

(section assembling); 

- way to mitigate the welding induced deformation (steel 
plates). 
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A computational simulation might aid researchers in 
planning changes in the shipyard without physical modifications, 
although other types of experiments will have to be performed 
so that to verify the actual results, as well. Because of this 
certain improvements demands practical trials to verify the 
actual results before the modification implemented. 

Looking forward to reduce the process timing, material 
waste, labour used and, consequently, waste of money, 
shipyards and design offices involve in wide-scale cooperation. 
For that reason, over the years researchers have implemented 
new ideas and carrying out various experiments in shipyards. 
This work is an example of that, showing a practical way to 
analyze the deformation in the ship panel after welding in the 
shipyard. Other examples of this kind of work can be found 
in [2], where a methodology allowing for the prediction of 
welding distortions and prevention of rework caused by, for 
instance, straightening, is presented. This way of thinking is 
also confirmed in the work [1], where experiments were used 
to understand the effects of welding on the structure, and in the 
work [7], which describes investigations of the weld geometry 
effects on the mechanical properties of undermatched welds. 

Generally, steel welding is the source of one of the most 
fundamental problems in the shipbuilding industry, as it leads to 
the deformation of steel parts due to residual stresses resulting 


from the interaction between welding parameters and the 
geometry of steel parts. Therefore this work aims at verifying 
a numerical way to predict the welding induced deformation 
after submerged arc welding during the normal production in 
the shipyard. 


WELDING PROCESSES IN THE 
SHIPBUILDING INDUSTRY 


According to [4], fusion welding can be described as 
a joining process which uses the base metal molten by high 
temperature to join elements.. 
The welding process has been used in many different types 
of industries due to: 
- the reliability of the structures assembled; 
- high speeds in mass production; 
- the possibility to use automated robots to perform the 
welding; 
- efficiency in constructing standardized parts through the 
use of automation; 
- high quality. 


Welding processes have disadvantages, as well, some of 
which are listed below: 

- complex quality control processes; 

- the need for skilled operators; 

- expensive and complex infrastructure when automated 
welding processes are implemented; 

- the possibility of substantial hazards and losses in the 
event a failure of the welded part during its life (the steel 
structures of civil buildings, vessels, oil rigs or motor 
vehicles). 


In order to avoid these disadvantages the welding technology 
needs special care about the process control. 

In shipbuilding, welding is extensively used because ofthe 
significant demand for steel or aluminium alloys to be used in 
structural elements, such as ship hulls, watertight bulkheads, 
decks, bottom, being typically all-welded thin plate structures, 
[3]. Therefore, during fabrication some deformations may occur 
that are caused by several factors, including the cutting of the 
plates and pieces, or the welding of the joints between them. 

There are several methods for connecting parts. Before 
commencing joining works, some issues must be considered: 
- the type of joint to be made; 

- the type of material to be welded; 

- the limitations of the production layout; 
- the thickness of the plates; 

- the purpose of the welding. 


Submerged Arc Welding (SAW) is widely used in the 
shipbuilding industry and, according to [7], it offers high 
productivity. It is a mechanized welding method that can make 
use of a single electrode or several ones at the same time. The 
electric arc burns underneath a layer of protective flux which 
melts in the vicinity of the arc and produces solid slag on the 
weld. The non-molten excessive flux 15 reused. 


A Few Characteristics 


To arrange the welding procedure and to get reasonable 
results this is necessary to understand a few characteristics 
and effects influencing the work pieces which will be welded. 
Some of them have been commented earlier, such as: the type 
ofthe welding process or workpiece geometry. However, other 
issues are equally important when the welding process is under 
evaluation. 


From this point of view, the following aspects should be 
treated as equally significant: 

Heat input 

Heat input can be described as one of the most important 
factors involved in welding and, according to [8], a significant 
variable related to the quality of welding, mainly with regard to 
the warping problem identified in welded products. Heat input 
affects the dimensions ofthe welded object, as it brings enough 
energy to melt the additional material and the workpiece. Such 
energy may generate local deformation, as indicated in [3]. Heat 
input is the relation between the welding energy used during 
the operation and the speed of the welding source along the 
process. This relation 15 measured in units of energy (Joules) 
per units of length (mm or cm). Together they are helpful in 
determining the amount of heat which will be applied to the 
metal in order to melt the base material (workpiece) and the 
fusible wire (additional material). 

Residual Stress 

The heat that 15 necessary to melt metal parts in order to 
join them reaches extremely high temperatures, inevitably 
causing certain physical phenomena inside the workpiece due 
to the difference between the molten pool temperature and the 
room one. As a consequence of this, some characteristics of the 
particular element may be affected due to the quick heating and 
cooling cycle that occurs during welding. 

Welding deformation 

Deformations generated during the welding process are 
strongly linked to the residual stress, and for many researchers 
they are an unavoidable problem present in every welding work. 
In the shipyards, according to [4] and [6], the most representative 
problem is the one of misalignment between blocks during the 
hull erection stage. These problems occur due to deformations 
caused by welding, which in turn cause inaccuracies in entire 
blocks and structural member misadjustment. These authors 
also claim that welding deformations do not only interfere with 
the aesthetical aspect, but also put the structural integrity ofthe 
ship at risk when in service. 

Hence, all these characteristics must be known during the 
welding process to minimize the problems once the ship is 
delivered. 


MODELLING 


The literature provides a multitude of methods developed 
over the years using the most diversified tools designed to 
identify, understand, correct or, at least, mitigate the undesirable 
effects that the aforementioned problem has on the workflow 
and other production aspects. Taking advantage of some of 
them may help eliminate certain inconveniences encountered 
in the workflow which, in turn, will improve the efficiency of 
the process. 

In order to counter these welding problems, FEM 1s usually 
used to simulate and investigate residual stress distribution 
along the geometry, and subsequently to suggest modifications 
to improve the process. 

For the same reason, statistical tools are used in the 
process analysis, especially to solve problems or even process 
management. This approach drives the research toward the 
understanding of the process, once a number of trials are 
necessary to be performed. In addition, this can lead to a better 
insight into the current situation, which may help minimize the 
unfavourable effects or increase the production capacity. 

The purpose of the experiments is to explain the process 
through its behaviours based on the actual results which are 
acquired from the process. The trials need to be carefully 
modelled, contemplating the major of the variables that take 
part in the analysis to know all effects of those variables in the 
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process. It means that all of the characteristics of the process 
should be evaluated in order to make sure that no inaccurate 
results are obtained. 

Itis very important to know which factors exist within the 
process that have a tangible effect on the final output. From 
the point of view of analysis, these are factors that contribute 
to the occurrence of the problem. To illustrate this, only a few 
of the various relevant factors are mentioned below regarding 
to the welding induced deformation: 

- welding energy (current and voltage); 

- welding speed; 

- welding type; 

- the thickness of the plates; 

- the boundary conditions during the welding process. 


In this case, experiments and evaluations were carried out 
during the normal production in the shipyard following the 
design of experiment (DoE) methodology. 


Defining the Problem 


According to [5], a number of experiments were carried 
out in Polish shipyard aiming at the prediction and mitigation 
of the welding-induced deformation which occurs after the 
butt-welding of 5 mm thick plates using SAW. The trials 
were developed on the stiffened panels with the following 
dimensions: width of 12 m and length of 10 or 12 m. These 
panels were used on the deck or in bottom area. These panels 
usually consist of 4 or 5 plates (depending on the length ofthe 
panel) and consequently 3 or 4 weld seams (12 m long). The 
plates were laid on the SAW work line to be joined by the weld 
seam, once their edges were already prepared. 

In this case, when the panels lay properly and flat on the 
welding line, the causes of such deformations were sought. 
Certain conclusions were drawn from the analyses and the 
shipyard staff came up with the idea to verify the welding 
gantry facilities used to perform butt welding on the first side 
of the panels. 

Fig. 1а provides a general view of the magnetic bed and 
the copper bar used for butt welding the first side of the panels. 
Fig. 1b is a simple diagram showing the forces (arrows) 
affecting the plates during welding. А Design of Experiment 
needed to be carried in order to figure out the effects of the 
welding parameters adding to the effect ofthe pressure induced 
by the copper bar and the magnetic bed. The first plate was 
positioned on the right side of the magnetic bed and then 
aligned along the edge and clamped by the magnetic blocks. 
The second plate was pulled towards the first edge plate until 
the span between them was around 0 to 5 mm, after this the 
second plate was also clamped by the magnetic blocks from the 
left side (Fig. 1b). Once properly positioned, the manual tack 
welding was performed by the operator, afterwards the operator 
puts down the tack welding using the pneumatic countersink. 
From this condition, the copper bar acted on the plate junction 
(1-ргоћ е) and then was possible to perform the first side butt 
welding using the SAW gantry. 

The copper bar has been used on SAW processes to 
guarantee that the molten material will not leak between the 
plates during the process, besides it works as ground for the 
electric procedure and holds the 2 plates at the correct position 
to perform the weld. 

The magnetic bed was used to fix the plates before welding 
and its working parameters were previously defined together 
with the experienced shipyard staff. The pressures were set 
up directly on the pressure board according to the values 
determined to be used during the trails (DoE — see Tab. 1.) 


76 POLISH MARITIME RESEARCH, No 2/2013 


= i Magnetics bed Ё 


ы 
L—tdGrT—3.8-3 ET! Lt ` QC—ÉEL— 
[NU | FA 


Magnetics bed 


Fig. 1. The welding gantry 


Data Processing 


The first idea for arranging the data was to evaluate the 
mean value for related points distributed along the entire 
panel. However, this plan was discarded due to the fact that 
each panel had 3 different weld seams connecting its plates. 
The deformation should be investigated by examining each 
weld seam with the intention of narrowing down the number of 
sources oferrors caused by external factors, and even excluding 
initial plate deformation. Thus, concentrating the measurement 
on areas near the weld seam could help to mitigate negative 
influence from these factors. 

Hence, a number of approaches were examined. On the one 
hand, we examined the mean value of deformation between 
points related to the weld seam and, on the other hand, the 
largest value of deformation on the way from the centre of the 
plate to the weld seam. 

Based on these examinations, two methods for data 
arrangement were analysed, as explained below: 

o Option 1 — Mean value arrangement 


This option is illustrated in Fig. 2 showing the diagram 
adopted for the purposes of obtaining the value. To facilitate 
the analysis, this arrangement will hereinafter be referred to 
simply as “Option 1”. 

This diagram consists of the mean value between two 
points (in the middle of the plate) neighbouring on a point on 
the weld seam. 
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Fig. 2. Mean value arrangement: Option 1 


Difference between 


4 7 “©. point 10 and 7 


Difference between | 
point 4 and7 


7 


Welding direction 


Fig. 3. Largest value arrangement: Option 2 


the value which determines the position of the 
deformation along the weld seam. This number 
varies from 1 to 3; 

the value measured in mm at point x; 

the index of point P. This number varies from 4 
to 24; 

the deformation calculated at weld seam j at 
point 1; 

the largest value among the three Do, у , values 
calculated at weld seam j. In other words, it is 1 
out of 3 values of deformation analysed during 
the data processing; 


Adopting this point of view, it is possible to realize that: 

- each weld seam has three different Do, ;, values (one for 
each edge and another one for the centre of the plate; 

- thelargest value among the three values measured is considered 
as an amplitude on the current weld seam Dop yea 

- consequently, each panel yields three values for data 
processing. 


* Option 2 - the largest value 


Fig. 3 illustrates the arrangement that takes into account the 
largest value of deformation around each weld seam. Similarly 
to Option 1, this arrangement will hereinafter be referred to as 
“Option 2”, simply to make the analysis easier. 

This diagram makes use ofthe largest value of the difference 
between individual points neighbouring on the weld seam. 
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Where: 
j — characterizes the weld seam. This number varies 
from | to 3; 
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1 — the value which determines the position of the Adopting this methodology, it is also possible to realize 
deformation along the weld seam. This number that: 


varies from 1 to 3; - each weld seam also has three different values (one for each 

P. — the value measured in mm at point x; edge and another one for the centre of the plate; 

x — the index of point P. This number varies from - the average value among the three values measured 15 
4 to 24; considered as an amplitude on the weld seam; 

Dop 5, — the deformation calculated at weld seam j at - consequently, each panel yields three values for data 
point i; processing. 

Doy веат) — the largest value among the three Do, 2, , values 
calculated at weld seam j. In other words, it is 1 Once the arrangement is defined, the values for plate 
out of 3 values of deformation analysed during deformation used during the data processing are showed in 
the data processing; Tab. 1. 


Tab. 1. Deformation values from two data arrangements 


Copper Bar Deformations 


(CB) agnetics (M) | Weld Seam Option 1 Option 2 
[bar] 
[bar] [mm] [mm] 


Heat Input (Q) 
[kJ/cm] 


шо | бој = | 6 | бо рк | бо | бо рк | бо | до | — | шә | до | бо | бо | — шә бо рк | бо | до | WIDOT — | бо | — орке | Оз | бо рк | Оз | бо | — | о | бо | Оз |t2 рә | — | 02 | F2] — 


78 POLISH MARITIME RESEARCH, No 2/2013 


RESULTS AND DISCUSSION 


Aided by the statistical software to compile the results 
obtained from the experiment, the values used generated two 
numerical equations. These equations make it possible to 
verify the results and the accuracy of each value provided by 
the models, by comparing it to the actual value of deformation 
measured on the steel ship panel in the shipyard. АП the 
variables were considered by the software, including the data 
arrangement (see Fig. 2 and Fig. 3), to get these 2 final equations 
and a number of consideration were taken into account to choose 
the models. These considerations are described below. 

For boundary conditions, the values for copper bar pressure, 
magnetic pressure and heat input (according to Table 1) were 
considered as well as the data arrangement to build the DoE 
(according to Figure 2 and Figure 3). The values from the plate 
edge were not considered and only real data (measured along 
the plate, according to the Figure 2 and 3) were inputted to get 
the equation to predict the welding deformation. 

The numerical equations are presented below: 
Option 1 — Mean value arrangement 


Do, = 185.70-13.49- Q-90.69:M +7.74-0-М (5) 
Option 2 – Largest value arrangement 


Пор,= 422.45 -35.87 - Q- 574.08. CB + 


о 


о 


(6) 
—113.80 -M + 52.41- Q- CB410.01- Q- M 


The values obtained from Option 1 demonstrate a standard 
deviation of approx. 6.1 mm, while the values obtained from 
Option 2 demonstrate a standard deviation of approx. 4.7. 
This difference exists because the data arrangements are 
different from each other. Therefore, an obvious consequence 
of that 1s that Option 1 presents a larger range of differences 
in deformation values than the one provided by Option 2 in 
trials illustrated in Tab. 2: 


Tab. 2. A comparison between the results from Option 1 and Option 2 


Amplitude of the 
differences in 
deformation values 
[mm] 


Limits 


Maximum 


Minimum 


Maximum 


Minimum 


A graphical comparison of these values is better at 
visualising the difference between the models. 

Fig. 4 (a) and (b) show that although the points are not 
fully distributed on the 45? straight line (which would be the 
ideal condition for the values predicted), they are concordant 
to a certain extent with the actual values. It is also obvious that 
the manner in which the points are laid down in Option 2 is 


Predicted vs. Actual 


30.00 
Actual 


Predicted vs. Actual 


40.00 50.00 


a) 
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Option 1 50.00 = 
Color points by value of 
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Fig. 4. A comparison between the actual and the predicted deformation values: a) Option 1 and b) Option 2 
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superior to the shape created in Option 1. This particularly refers 
to the points located at the extremities ofthe curve (the larger 
and smaller values obtained in the arrangement). In Option 2, 
they are spread along the line, whereas the points in Option 1 
are concentrated near the mean values. 

Another method for verifying the predicted deformation 
values for both the cases studied is shown in Fig. 5. 

The cube chart makes it possible to verify all the predicted 
deformation values across the entire spectrum ofthe variability 
of the main factors. Highlighted in Fig. 5 (a) and (b) are 
the maximum (the black rectangle) and minimum (the grey 
rectangle) results for the deformation predicted for both the 
models. The best condition for work (i.e. the condition which 
renders the lowest deformation values according to those 
models) is with the following values: heat input = 10.44 kJ/cm, 
copper bar pressure — 0.45 bar, magnetic pressure — 2 bar. 

According to these results, Option 2 demonstrates better 
results than Option 1 due to the smaller amplitude of the 
residual variation. 

After the results of these two options were verified from the 
statistical and dimensional points of view, it was confirmed that 
the model described by Option 2 yielded better outcome than 
the one described by Option 1. There are a number of reasons 


a) 
Design-Expert® Software 
Factor Coding: Actual 
Option 1 

X1 = A Heat Input 
X2 = C: Magnetics 

ХЗ = B: Copper Bar 


that justify these results, one of them being leaving out the 
copper bar effects of the model. 

Moreover, another aspect of the approach adopted in this 
paper is the fact that the models developed through the DoE 
account for the average of the input values, 1.е. the models 
operate based on the average deformation between the weld 
seam values. If we compare the predicted values to the average 
values of the weld seam the standard deviation becomes smaller, 
which in turn means that the model is capable of properly 
reflecting most of the results. 

Therefore, for future evaluations using this methodology, 
it would be reasonable to evaluate the welding deformation of 
5 mm thick steel plates. 


SUMMARY 


After the work done in the shipyard, the numerous trials 
performed on the panel welding line, the measurements taken 
and the result analyses carried out, a number of interesting 
conclusions can be drawn from the investigation, which can 
be formulated as follows: 

1. Unlike the main methodology making use FEM simulation 
as proposed by several authors, the welding deformation in 


Cube 
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Fig. 5. Predicted deformation values: a) Option 1 and b) Option 2 
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butt joints on 5 mm thick steel plates can also be predicted 
by a statistical method, as shown in this paper. 


. An analysis subject to statistical principles may yield 


different results depending on how the data are inputted, and 
this characteristic became an important issue encountered 
in the course of this paper. Because of this, specific data 
arrangements were chosen and their results verified with 
the intention to derive a numerical equation capable of 
predicting, with the highest accuracy possible, the values of 
deformation similar to the ones measured in the shipyard. 
Surprisingly, from each data arrangement a different 
equation could be derived and, sometimes, such an equation 
did not take into account the influence of the all of the 
factors present. 

The numerical equation obtained from the statistical 
software, through the data arrangement described by 
Option 2, rendered reasonable results, as compared to 
the data collected in the shipyard. Therefore, it can be 
considered a satisfactory method for predicting the welding- 
induced deformation in 5 mm thick steel panels. 


It is also worth noting that all kinds of experiments, 


especially those performed under strict conditions dependent on 
time and equipment availability, are subject to the influence of 
the existing sources of errors. Such sources should be identified 
and, whenever possible, avoided. For future projects related to 
analyses through experiments, the following recommendations 
should be taken into account: 


l. 


Measurements: the same person, the same equipment 
and the same conditions should be used to perform the 
measurements. 

Feedback: throughout the task, there is a need for intensive 
communication with the workers responsible for the process 
to make sure they have sufficient understanding of the 
experiment. These workers will be in charge of the quality 
control; therefore, if during the trials something goes wrong 
they will be able to introduce the required corrections and 
fix the problems as they occur. 

Feedback: similarly, the key person (tasked with rendering 
assistance whenever necessary) must be identified before 
the start of the experiment. 

Noise in the model: verification of the initial condition 
(deformation) ofthe plates before the commencement ofthe 
welding process. It can help the researcher to sort the noise 
out of data processing, improving the overall results. 


A few of the aforementioned errors may have some 


minor influence on the final results presented. However, the 
methodology used was developed taking as careful measures 
as possible to ensure its reliability. 


A number of extra trials should be performed in order to 


verify the actual accuracy of the model. 
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ABSTRACT 


e 
©, 


The article discusses a problem of performing a mission and controlling a team of unmanned 
underwater vehicles which perform a task of water region search. The team of underwater 
vehicles has been modelled as a multi-agent system, which made it possible to use the 
mechanism of multi-agent systems for controlling the team of cooperating underwater 
vehicles, with each vehicle being treated as an individual agent. Solutions to cooperation 
problems in multi-agent systems were proposed, including the presentation of the system 
structure and the data exchange method. Issues relating to cooperation actions in the team 
of underwater vehicles were described and solved using the theory of games and the vector 


quality index. The process of negotiations between the vehicles was modelled as the multiplayer cooperative 
game, which secured successful realisation of the task performed by the team of vehicles with respect to time, 
which represents the economic aspect of the task. The optimal solution was selected based on the vector 
quality index, which provided opportunities for minimisation of costs having the form of energy used by the 
entire team, and maximisation of profits in the form of time needed for completing the water region search 
task. The proposed method for dividing the action space in the water region search tasks performed by the 
team of underwater vehicles consists in assigning to each of cooperating underwater vehicles a task the 
realisation of which contributes to the realisation of the globally defined task. Finally, the results of water 
region search tasks performed by a team of cooperating underwater vehicles are presented and discussed 
taking into account possible vehicle failures happening during task realisation. 


Key words: unmanned underwater vehicles; multi-agent systems; control; multiplayer cooperative games 


INTRODUCTION 


Nowadays, research activities concerning generally 
understood robotics are frequently oriented on the use of 
teams of cooperating robots performing tasks which cannot 
be performed, or will be performed less efficiently by a single 
robot. Such a situation can be observed in underwater robotics 
where concepts are developed for certain tasks to be performed 
by a team of underwater vehicles. These tasks most often 
include exploration and inspection activities, which in general 
can be done by a single robot but a team of robots can do them 
faster and more effectively. Another kind of tasks of this type 
are transport activities in which overall dimensions and the 
shape of the transported structure make the task to be only 
done by a team of cooperating vehicles. 

A team of cooperating robots reveals numerous advantages, 
such as high resistance to failures, for instance, which makes 
it possible to perform the task even if one of robots beaks 
down. Moreover, in some cases the decomposition of the task 
to be performed by the team of robots allows us to use robots 
with simpler structure than that of a single robot performing 
the same task. However, the problem of controlling the team 
of cooperating robots is much more difficult than controlling 
a single vehicle. In case of a team of robots, additional issues 
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are to be discussed, including communication between the 
vehicles, coordination of actions, performing negotiations, 
and many others. 

The computer sciences make more and more frequent use of 
multi-agent systems, especially to solve problems of distributed 
nature, or those requiring complicated calculations. The multi- 
agent system comprises a number of agents which communicate 
and cooperate with each other to perform a common task. The 
mechanism ofthe multi-agent system can be easily adopted to 
control the team of cooperating underwater vehicles in which 
each vehicle is treated as an individual agent. 

Below presented is a concept of a multi-agent system 
which was used to perform a task of water region search 
by a team of cooperating underwater vehicles. The team is 
assumed to consist of homogeneous and cable-free underwater 
vehicles equipped with instruments for navigation, technical 
observation and communication used for data exchange 
between the vehicles. The underwater vehicles have docking 
stations situated at one point of the water region, where they 
can load their internal energy sources, transmit the collected 
data and receive commands to start work. The shape of the 
searched water region 1s defined by the vertexes of the polygon 
enclosing the region. The water region is divided by the superior 
system into areas, each of which is to be searched by a single 


underwater vehicle. The information on the area to be searched 
and the corresponding passing trajectory is obtained by each 
vehicle from the docking station, together with the order to 
start work. 


AGENT 


The literature on the subject does not deliver any precise 
definition of an agent. Generally we can say that the agent is 
an element which is used for creating decentralised systems 
[2]. Another definition says that the agent is an object which 
can both observe the environment and act on it [15]. The most 
precise definition describes the agent as a system situated in 
a given environment which can make use of certain resources 
and be able to act autonomously in order to achieve assumed 
goals [4, 13, 14]. To make the underwater vehicle fulfil the 
above definition, it has to be equipped with sensors to record 
input signals coming from the environment, and effectors 
using which it can act on the environment. Additionally, is 
should have sort of internal awareness which will secure its 
autonomy manifesting itself as making decisions on actions to 
be performed to achieve the assumed goal. We can say that the 
agent can adapt its behaviour to a current situation. It should 
be noted, however, that the goal to be performed by one agent 
can differ from that assigned to another agent, and in case of 
common task realisation the request coming from the other 
agent can be realised only when it is in line with the interest 
of the asked agent. 

To describe formally the action of the agent we introduce 
the following notations: S = {s,, s,, ...} represents possible 
environment states, and A = {a}, a, ...} possible actions of the 
agent. In this situation the action can be meant as the function 
[13, 4]: 


action: 5" > A (1) 


where: 


S* — the selected sequence of environment states. 


The environment can be defined as the function: 
environment: S x А — S (2) 


The above notations mean that the undertaken action results 
from the environment state, while the environment modifies 
its state due to the action of the agent. 

We should also define a function which allows the agent to 
read input signals. This action can be written as the function: 


perception: S — P (3) 


which maps the state of environment into perception P. 
Consequently, the action function takes the following form: 


action: P* > A (4) 


where: 
Р* — the perception sequence. 

Assigning the property of autonomy to the agent consists 
in defining the function of actions which are to be initiated by 
the agent depending on input signals composing the perception. 
To do so, we should select a method of implementation of the 
decision function, at the same time selecting an appropriate type 
ofagent. In case ofunderwater vehicles it seems a best solution 
to use the hierarchical layer architecture, the advantages of 
which include no need to define the decision making mechanism 
and simplicity of action. A popular sample of this type of 
architecture is the InterRRaP architecture (Fig. 1) [11, 12]. It 
consist of a series of layers, the lowest of which is responsible 
for reflex actions, the middle layer for planning, and the top layer 
for cooperation with other agents [4, 11, 12]. Each layer has its 
database which represents the state of environment attributed 
to this layer. The flow of control in this architecture is realised 
in two directions: first the control signal is transmitted from 
lower to higher layers, and then from higher to lower layers. If 
the lower layer, making use of its database, is not able to take 
a relevant action, it sends the signal to the higher layer. That 
means that the higher the layer, the more rarely it is activated. 


ENVIRONMENT 


Input data from sensors 


Reflex action layer 


Planning layer 


Action undertaken 


— 


Cooperation layer 


Social knowledge 


Fig. I. Flow of control in InteRRaP type architecture [4] 
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According to the adopted assumptions, a single underwater 
vehicle is an agent whose task 15 to search the water region in 
cooperation with other vehicles. Therefore the lowest layer of 
the agent is responsible for avoiding obstacles [17]. The middle 
layer is responsible for precise motion ofthe underwater vehicle 
along the given trajectory [16], while the role of the highest 
layer is to cooperate with the remaining vehicles in the team. 


MULTI-AGENT SYSTEM 
Structure of the multi-agent system 


A team of cooperating underwater vehicles in which each 
vehicle is an agent composes a multi-agent system. The agents 
in this system communicate with each other, at the same 
time preserving the autonomy of their operation and decision 
making, which is reflected in actions undertaken by particular 
agents. Such an interaction between the agents leads to the 
modification of knowledge and actions performed by individual 
agents based on the observation of the behaviour of other agents 
and information transmitted in the communication process. 
An individual agent is equipped with the strategy leading it to 
achieving the assumed goal. However, being part of the multi- 
agent system it has to put in practice the strategy of cooperation 
with other agents, especially when this strategy facilitates 
reaching the goal common for the entire team. Cooperation 
and competition between the agents in a multi-agent system 
is its key feature [9]. To make it possible, we should define 
the structure and formulate protocols of interaction to provide 
opportunities for describing message exchange scenarios and 
allow making decisions. 

The multi-agent system can have a hierarchic structure, 
in which a decision about future actions is made by the team 
leader, or a flat structure, in which all agents are equal and 
decisions are made as a result of negotiations. The advantage 
of the hierarchical structure is high speed of decision making, 
as the time-consuming process of negotiations can be avoided 
here. Its disadvantage 15 relatively poor use ofall team potential, 
asthe team leader can make decisions based on the information 
on the environment and other agents which is insufficient to 
make the right decision. The problem also arises when the 
current leader breaks down and mechanisms of selection of 
a new leader are to be applied. In case of flat structure this 
problem does not exist and decisions are made collectively. 
Unfortunately, decision making in this structure is a very time 
consuming process. 

Either ofthese two structures can be successfully used in the 
task of water region search realised by a team of autonomous 
underwater vehicles. However, since the water environment 
is very dangerous and unfriendly to technical equipment, we 
should assume that failures of vehicles can happed relatively 
often. Moreover, search tasks can be performed in relatively 
long times, which means no limits concerning the speed of 
decision making processes. That is why the flat structure seems 
more appropriate and better promises success in the realisation 
of tasks of the above type. 


Data exchange 


The essence ofthe multi-agent system 15 coherent behaviour 
of all agents composing the team. Difficulties in the realisation 
of this task result from the fact that multi-agent systems have по 
central control, and securing coherent action of particular agents 
requires data exchange via mutual communication. Therefore 
we should assume that each agent has ability to send, receive 
an analyse the data, and that there is a well-defined protocol 
of interaction, the same for all agents. 
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The interaction protocols are divided, in general, into 
coordination protocols and communication protocols. The 
coordination protocols are used for periodical exchange of 
information on current condition of particular agents, including 
redundancy avoiding, while the cooperation protocols are used 
for settling cooperation between particular agents, i.e. task 
decomposition and distribution. The action of both protocols 
requires a relevant mechanism of data exchange. An interesting 
solution in this context seems to be the interaction mechanism of 
Blackboard type, the idea of which consists in the existence of 
common memory, frequently bearing the name of a blackboard, 
in which agents can write and read data. This way each agent 
has access to the data recorded by other agents and can, at an 
arbitrary time, write and read data in the blackboard. The data 
stored in the blackboard are valid for certain time, or until they 
are updated by more current data. Since the adopted structu re 
is flat without a defined leader which would be responsible for 
data storage, the common memory will be stored by each agent. 
In other words, each agent will have a copy of the information 
stored in the blackboard. In this case data modification in the 
blackboard will consist in transmitting new data by one agent 
to all remaining agents, after which the data will be updated 
in the blackboard. 

Since the underwater vehicles perform their tasks in the 
water environment, the common transmission medium of all 
agents is water, and the transmitted messages will have the form 
of mechanical (sound) waves. In this context it seems reasonable 
to use a token ring mechanism, which was successfully used in 
computer networks. The action of this data exchange mechanism 
makes use of the exchange of authorisation, or a so-called token. 
The vehicles are numbered in accordance with the increased 
distance of the area to be searched by them from the docking 
place. The agents, one by one, receive a token which allows 
them to pass a message, and then pass the token over to the next 
vehicle. The last vehicle passes the token to the first one. The 
agent which wants to pass the message waits for a free token, 
and having received it changes it into the occupied token and 
passes the data block bearing the name of a frame. 


Cooperative action 


The cooperative action in multi-agent systems can be 
described as a multiplayer cooperative game. The cooperative 
game will be an ordered pair having the form [1]: 


Г = (№, v) (5) 
where: 
N= {1,..., n, ..., N} - the set of player numbers 
V:2N Ri! — the characteristic function of the game. 


The characteristic function of the game assigns to each set 
S CN, where S represents all possible coalitions in the team 
of agents including those composed of single agents, a total 
maximal secured payment v(S) which can be obtained by the 
members of coalition S irrespective of the action undertaken 
by the remaining players. It is additionally assumed that 
0(@) = 0 and v(-) is super-additive i.e.: 


A(S) + 0(Т) xwS UT); SS TC N; ба T - 2 (6 
A solution to the multiplayer cooperative game can be 


presented as the solution to the multi-criteria optimisation of 
the following type: 


(X, F, R) (7) 
where: 
X — the nonempty set of feasible solutions, 
F:X—N - the vector quality index, and R defines the relation 
of domination. 


When analysing the cooperation between the agents 
performing the task of water region search, the set of feasible 
solutions is the set of all possible coalitions which can be 
formed between these agents. The criterion function assigns to 
each feasible solution x € X its numerical valuation, important 
for selecting the best solution. In search-type tasks we tend 
to minimise the time of task realisation, which 1s directly 
connected with the minimisation of the length of the path to 
be covered by each vehicle. Since the cooperation game 15 
played when one of agents cannot do the assigned task and 
another agent, or coalition of agents, is to take it over, we will 
want the profit from the task done by this vehicle or coalition 
of vehicles to be the maximal, the cost borne by each vehicle 
to perform the supplementary work to be the minimal, and the 
coalition performing the supplementary work to be the largest. 
The above assumptions would allow the team of vehicles to 
do the water region search task at the shortest possible time. 
Consequently, the quality index can be defined as: 


F(x) = {Z, K, L} (8) 


where: 

the profit which can be obtained as a result of 
supplementary action undertaken by the coalition of 
agents, 


K — the cost borne by each agent composing the coalition 
formed to do the supplementary work, 
L — the number of coalition members. 


Let W, represent the cost which a given agent has to 
bear to reach the place from which the supplementary action 
begins. This cost can be identified with the length of the path 
between the last point of its own trajectory and the point from 
which the supplementary action begins. Additionally, let W} 
represent the cost connected with covering the remaining part 
of the trajectory within the framework of the realisation of the 
supplementary task, which can be identified with the length of 
the trajectory to be covered in this supplementary action. In this 
case the cost K „a is the maximal cost borne by a single vehicle 
performing the supplementary task with the longest approach 
path, which can be written as: 


К max= Wr + max (Wp ;) (9) 
ieN 
where: 
i = the agent number from within a set of all efficient 
agents. 


When a coalition between a number of agents is formed, 
the task will be performed simultaneously by all coalition 
members. Therefore the cost borne by a single agent composing 
the coalition will be the sum ofthe cost borne for realisation of 
the remaining trajectory part and the cost of approach, taking 
into account the number of agents composing the coalition, 
which can be written as: 


= ies (10) 


The profit Z obtained by the given coalition of agents will 
be defined as: 


Z-K,.-K (1) 


which means that it 15 the difference between the maximal cost 
of realisation of the supplementary work by a single agent with 
the longest approach path and the cost to be borne by the agent 
as a coalition member. 


max 


A set of all values of the criterion function F which it takes 
over a set X creates the criteria space Y: 


Y= {у= Ех) є 9»|xe Xi 


being the map of the set X in Jt. 

To select a desired solution x from the set X we have to 
compare different solutions. This comparison is realised by 
comparing numerical values of particular estimations F,(x) 
attributed to the compared solutions. Therefore we can say 
that the decision whether solution x, is better than solution x, 
is made by comparing their numerical images F(x,) and F(x,). 
We should define the relation R С Y x Y, i.e. the relation of 
domination which will make it possible to decide whether the 
estimation ofthe image F(x,) = y, is better than that of F(x,) = у, 
or not. 

The relation R can be defined using the propositional 
function in the following way: 


Ко= (у у) е Ух У (у, у) =} 113) 


where the propositional function q: Y x Y — {0,1} has the 
form: 


(12) 


Ф(у,, Y2) = y, is better than y, (14) 


In the analysed case, for the quality index having the form 
of (3.16) the relation of domination can be defined in the 
following way: 


Ф(у„у) =L, > L; (15) 
Ф(у, У) = 7, 2 Z; 
and (16) 


Z2 IL] max Z({n}) л Z, 2 IL, | max Z({n}) 
nes, neS, 


Py У) = K, SK, 
and 


Kui — Wi pi 22Wip; ^ КУ m" №, pi > 2W,p; 


(17) 


Condition (15) says that the coalition performing the 
supplementary task is to be as large as possible. Condition 
(16) says that the profit to be obtained as a result of coalition 
formation 15 to be the maximal, and at least L times as big as the 
maximal profit of all coalition members, where L is the number 
of coalition members. Condition (17) says that each coalition 
member is to have an attributed search trajectory with the length 
larger or at least equal to two lengths of the approach path to 
the point at which the supplementary action starts. 


DIVISION OF ACTION SPACE 


When a team of cooperating underwater vehicles performs 
a task, a problem which is to be decided upon 15 how to divide 
the action space. The division of the action space consists in 
allocating an individual task to each of cooperating underwater 
vehicles is such a way that all individual tasks compose the 
globally defined task. In case of water region search task 
performed by a number of underwater vehicles, the division 
of the action space will consist in physical separation of water 
region areas to be searched by a single underwater vehicle. 
Searching all areas of the water region divided in the above 
way will be equivalent with the realisation of the search of the 
entire water region. 


POLISH MARITIME RESEARCH, No 2/2013 85 


A solution to this problem consists in geometrical division 
ofthe water region treated as a polygon into smaller polygons, 
being separate areas for which probing lines and underwater 
vehicle passing trajectories are then determined [3, 7, 8, 10]. 
However, this solution has some disadvantages. For instance, 
in this solution the time needed for reaching the searched area 
by the underwater vehicle, a parameter of certain importance, 
cannot be estimated. Moreover, when determining probing lines 
and passing trajectories for each of the separated areas, it may 
happen а ће border between the adjoining areas that some area 
parts are omitted, or are searched twice or more times. These 
disadvantages can be eliminated by introducing another method 
of action space division, which consists in determining first the 
probing lines and passing trajectories for the entire searched 
water region, and then dividing the trajectories into segments, 
each of which to be covered by a single underwater vehicle. 
Such a solution has certain advantages, the most important 
of which is elimination of time-consuming procedures of 
geometrical division of the water region and substituting it by 
a simple iterative algorithm of trajectory division into segments. 
Moreover, after combining together, the data recorded by all 
underwater vehicles performing the water region search task 
have the same form as those obtained, hypothetically, by 
a single underwater vehicle. The proposed solution also secures 
uniform distribution of work over all underwater vehicles and 
provides opportunities for covering the entire area of the water 
region by technical observation means without multiplying the 
observation. Figure 2 shows differences between dividing the 
action space on the geometrical basis and using the passing 
trajectory division. 


a) 


b) 


Fig. 2. Division of action space: a) geometrical division of water region, 
b) division of passing trajectory. 


The basic idea of the proposed solution consists in 


dividing the full trajectory worked out for the entire searched 
water region into segment explored by particular underwater 
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vehicles. In order to meet the condition that all underwater 
vehicles taking part in the search are similarly loaded the 
trajectory 15 to be divided in such a way that each segment has 
approximately the same length, including the paths of vehicle 
approach to the starting point of search operation and return 
to the docking place. 

The task consisting in trajectory division into segments can 
be defined in the following way: letn > 2 represent the number 
of segments into which the trajectory SK is to be divided. 
We assume the existence of such points А, belonging to the 
trajectory SK that: 


БА]+ X [Aa 
i-l..n-2 


+|Аз-к|=БК| an 
and that 

БА |+, = ЈАЈА |+4 = 
(19) 


=... = АЈА 


+d; =... = 


А.К +da 


where: d; is the total length of the shortest path of approach to 
point A, and return from point A, ,, 1.e.: 


d; - A, 


1+1 


4 А, 9 (20) 


1+1 


Condition (18) secures that the entire trajectory will be used 
when dividing the trajectory into segments, while condition (19) 
means that each trajectory segment, together with the approach 
and return paths, will have the same length. 


RESULTS OF INVESTIGATIONS 


In previous sections a solutions was proposed which 
allows controlling a team of underwater vehicles cooperating 
within the framework of the multi-agent system. The control 
algorithms making use of the multi-agent system and combined 
with the algorithm of control along a given trajectory 516 were 
implemented on the computer, which provided opportunities 
for performing their tests. 

As а test sample, the task of Gdynia harbour search was 
performed by a team consisting of eight cooperating underwater 
vehicles. It was assumed in the test that the underwater vehicle 
2 broke down after covering 3 947.6 [m] of the trajectory 
attributed to it. According to the adopted solution the remaining 
underwater vehicles started negotiations to settle which vehicle, 
or coalition of vehicles will complete the task started by vehicle 
2. After negotiations, the vehicles presented costs of realisation 
of the supplementary task, see table 1. 

As a result of the performed negotiations, profits from 
particular coalitions were calculated, and then coalitions 
securing highest profits were selected in groups of different 
numbers of coalition members (table 2). It was checked for 
each coalition whether the conditions (15) (16) and (17) are 
met. Based on that it was concluded that the most profitable 
solution will take place when the supplementary task is taken 
over by a coalition of agents consisting of vehicle 3 and vehicle 
5. Since the coalition consists of tvo members the profit will be 
2.96 times as high as in case of single vehicle coalition, while 
the cost of search will be 5.75 times as high, compared to the 
cost of approach of vehicle 3, and 3.76 times as high as for 
vehicle 5. In case of larger coalitions, condition (17) is not met 
for one or more vehicles. For coalitions composed of more than 
3 vehicles condition (16) is not met either. Table 3 collects the 
numerical data describing the performed task, while Figure 3 
shows the passing trajectories of each underwater vehicle 
against the electronic map of the searched water region. 


Tab. 1. Presented costs of supplementary action to be undertaken afier failure of vehicle 2 


Vehicle number 1 3 


4 5 


6 7 


8 


Presented cost 1 100.6 526.4 


1 042.9 769.2 1 404.1 1417.3 1942.7 


Tab. 2. Data settled in the negotiation process after vehicle 2 break-down 


Coalition 


3; 


1416.3 


Cost borne by a single vehicle 
6336.9 


B, 5} 


4200.15 


3553.0 


{3, 4, 5} 


5036.9 


2716.3 


{1, 3, 4, 5} 


5440.8 


2312.4 


{1, 3, 4, 5, 6} 


5622.5 


2130.74 


{1, 3, 4, 5, 6, 7} 


5741.4 


2011.8 


{1, 3, 4, 5, 6, 7, 8} 


5751.2 


2002.0 


Tab. 3. Data of the performed task of Gdynia harbour search 


Harbour area [012]: 1714708.2 


Number of vehicles in the team [-]: 


8 


Shoreline length [m]: 15334.7 


Vehicle 1: 


path length [m]: 


9 676.2 


realisation time [min]: 


80.65 


Water region depth [m]: 


Vehicle 2: 


path length [m]: 


3 947.6 


realisation time [min]: 


33.0 


Distance from bottom when 
performing the task [m]: 


Observation zone width [m]: 


Vehicle 3: 


Vehicle 4: 


path length [m]: 
realisation time [min]: 
path length [m]: 


12 640.8 


9 867.1 


realisation time [min]: 


82.5 


Distance between probing 
lines [m]: 


Vehicle 5: 


path length [m]: 


13 044.1 


realisation time [min]: 


108.7 


Course of probing lines: 


Vehicle 6: 


path length [m]: 


9 437.1 


realisation time [min]: 


78.7 


Total path length [m]: 42494.4 


Vehicle 7: 


path length [m]: 


9 647.7 


realisation time [min]: 


80.4 


Total time of task realisation 


by a single vehicle [min]: 354.12 


500 1000 


Vehicle 8: 


1500 2000 
Fig. 3. Passing trajectories realised by each of underwater vehicles composing the team which performed the task of Gdynia harbour search 


path length [m]: 


9 879.3 


realisation time [min]: 


2500 


82.5 
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Based on the performed tests we can conclude that the use of 
many vehicles for searching a water region leads to considerable 
reduction of time needed for doing this task. It is noteworthy, 
however, that this time does not decrease proportionally to the 
number of the used vehicles, which results from the fact that 
the time of approach of each vehicle to the area of its operation 
is to be taken into account. Assuming no vehicle break-down, 
a single vehicle needs 129.2 [min] for covering the trajectory 
having the length of 15503.9 [m], while the cooperation of 
five underwater vehicles extends the trajectory by 7623.6 [m], 
which is the total approaching distance of the vehicles to the 
area of operation. However, at the same time the total searching 
time is 3.3 time shorter and is reduced to 39.4 [min]. Assuming 
that one vehicle breaks down during task realisation, a single 
vehicle needs 354.12 [min] for covering the trajectory having 
the length of 42494.4 [m], while the team of eight cooperating 
vehicles extends this distance nearly twice, to 78139.9 [m], but 
the time of task realisation is equal to 108.7 [min], i.e. 3.3 times 
faster than it would be done by a single underwater vehicle. 
Another advantage of coalition which should be mentioned 
here is that the search task for the entire water region can be 
successfully completed, which would not be possible in case 
of one vehicle performing the search and its subsequent failure 
during task realisation. 


CONCLUSIONS 


Controlling a team of autonomous underwater vehicles 
is a complicated issue, especially in case when the motion of 
vehicles is subject to limitations oriented on minimisation of 
the task realisation time. The problem becomes even more 
difficult when the vehicles have to cooperate to perform 
supplementary tasks after a failure of one or more vehicles 
taking part in the search action. Taking into account water as 
unfriendly environment and difficulties connected with precise 
realisation of the set trajectory, as well as with navigation 
and communication, we should conclude that the problem 
of controlling the team of underwater vehicles is extremely 
nontrivial. Its solution needs an untypical combined approach, 
which will link the knowledge from within the fields of 
automation, computer science, communication, and sociology. 
All these features can be found in multi-agent systems which 
are still in the phase of tests an definition of principles. Their 
basic idea secures cooperation of a team of robots acting 
to achieve acommon goal which is the realisation of the 
formulated task. 

The performed tests have revealed that controlling the 
team of underwater vehicles with the aid of a multi-agent 
system provides opportunities for efficient realisation of the 
task defined as searching a given water region. Adopting the 
above concept to the underwater search task resulted in doing 
the task in much shorter time than in case of a single vehicle 
performing the same task, simultaneously securing higher 
chances of success of the mission even when one or more 
underwater vehicles break down. The proposed solutions 
permit more efficient distribution of individual tasks, being 
parts of the global task, over individual vehicles composing 
the team. The proposed communication solutions make it 
possible to exchange efficiently the information between the 
agents, which secures fluent realisation of both coordination 
and cooperation actions. Moreover, defining the cooperative 
actions as the multiplayer cooperative game with vector criteria 
function allows us to use multi-criteria optimisation to solve 
problems of supplementary actions to be undertaken after 
failure of one or more agents. The adopted quality indices and 
the defined functions of domination secure the realisation of the 
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water region search task by a team of underwater vehicles in 
minimal time and at minimal costs. The proposed algorithms of 
agent’s operation along with the defined actions of cooperation 
performed during the negotiations secure proper behaviour of 
individual agents and, consequently, of the entire team. 

It is noteworthy that the combined use of a number of 
cooperating underwater vehicles in water region search tasks 
result in shortening the time needed for doing the task. However, 
the time reduction is not proportional to the increasing number 
of the used vehicles, which is connected with extra time 
expenses to be covered by each vehicle to reach the area of its 
action. Taking this into account when dividing the action space 
makes it possible to assign optimally partial tasks, while at the 
stage of cooperation actions it allows to choose the coalition 
of vehicles performing the supplementary action in sucha way 
that both the time of operation of the entire team and cost of the 
supplementary operation are minimal. The adopted solutions 
make it possible to determine the minimal time needed for 
performing the search action by the team of underwater vehicles 
in the given water region, even after assuming the failure of one 
or more vehicles before the start of the search. Moreover the 
proposed solutions make it possible to determine the optimal 
number of underwater vehicles which will perform the action 
in the given water region, thus minimising both the time of 
search and the number of vehicles. 
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ABSTRACT 


Transport of mineral concentrates is associated with high risk of liquefaction of the cargoes. 
This type of goods containing water may liquefy and shift dangerously across the cargo 
compartments. Therefore, it is rational to limit the moisture content of the cargoes which 
may liquefy. To prevent sliding and shifting of ore concentrates in storage biodegradable 
materials composed of starch are added to the mineral concentrates. The influence of 
adding starch materials to the mineral concentrates on its parameters determining ability 
for safe shipment was assessed on the basis of determination of the following parameters: 


permeability, cohesion and internal friction angle of concentrates. 


Keywords: mineral concentrates; liquefaction; shifting; starch materials 


INTRODUCTION 


Bulk cargo shipping has been used for many years to reduce 
the cost of sea transport and the transport of bulk cargoes is 
a vital component of international trade. Such trade requires 
a sufficient volume of cargo suitable for bulk handling and 
hence justifies a tailored shipping operation. The five major 
dry bulk cargoes are: coal, mineral concentrates, grain, bauxite 
and rocky phosphate, and the trade in bulk increases each year 
[1]. 

Mineral concentrates and other similar fine — grained 
materials are transported by sea in a wet state. The process has 
become increasingly associated with high risk of liquefaction 
of the cargoes. The International Maritime Organization 
(IMO), recognising that some losses had occurred due to 
improper loading, issued codes of practice for such operations 
[2, 3, 4]. 

To each of mineral concentrates liable to liquefaction 15 
assigned a Transportable Moisture Limit (TML) representing 
the maximum moisture content of such cargo to be considered 
safe for carriage in ships. The cargoes prone to liquefaction 
should never be carried without checking their moisture 
content. The International Maritime Solid Bulk Cargo Code 
(IMSBC Code) lays down that a certificate stating the relevant 
characteristics of the material to be loaded, including also its 
TML value, should be provided at the loading port. The Code 
provides information how the moisture content of mineral 
concentrates can be tested and assessed. 

The liquefaction is caused by migration of moisture. The 
moisture migration occurs when vibration at sea progressively 
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makes the water contained in a cargo to migrate down to the 
bottom ofa hold. The process might eventually lead to forming 
a free surface triggering dangerous violent shifting the cargo 
from side to side. The resulting much wetter bottom layer may 
therefore be prone to liquefaction and provoke instability of 
the entire cargo even though the average moisture content of 
the whole cargo remains lower than the TML [5, 6]. 

In general, finer material will liquefy (due to rise of water 
pressure in pores) more easily than coarser grained cargo 
because the water moves more easily through the larger pores 
between coarser particles. As a result, the pore pressure in 
coarse material dissipates more rapidly or does not build 
up at all. This is a generalization because many cargoes аге 
not all fine or all coarse but usually consist of a mixture of 
particles of various size. They have a wide distribution of size 
of particles. 

Natural and synthetic polymers (and their respective 
composites) are currently extensively used as biomaterials 
and in biodegradable applications because they are available 
in a wide variety of composition and their properties can 
be tailored to meet specific demands. Differences in the 
manufacturing process and in the composition of the raw 
material can significantly alter the functionality of the final 
product [7, 8]. 

In recent decades a major position has been taken by 
the materials made of synthetic polymers processed from 
crude oil. The cost of their production is low while their 
functional properties are very good and so 15 their resistance 
to environmental influences. Their durability regarded as 
a favourable feature, has become a serious ecological problem 


[9]. As а result, attempts have been made to solve the problems 
by introducing biodegradability into polymers through slight 
modifications of their structure. 

Inthe last few years it has been noticed an increasing interest 
in biodegradable polymers derived from annually renewable 
feedstock such as starch and starch-based materials [10, 11, 
12]. Starch is a polymer which occurs widely in plants. In all 
plants starch is produced in the form of granules. In general, the 
linear polymer, amylose, makes up about 20% of the granule, 
and the branched polymer, amylopectin, the remainder. The 
starch molecule has two important functional group, the -OH 
group that is susceptible to substitution reaction and the C-O- 
C bond that is susceptible to chain breakage. The hydroxyl 
group of glucose has a nucleophilic character. By the reaction 
of its -OH group, modification of various properties can be 
obtained. Cross-linking or bridging the -OH groups changes 
the structure into a network thus increasing the viscosity and 
reducing water retention. 

Chemical and physical properties of starch have been widely 
investigated because it is easy to convert it into a thermoplastic 
and then to use it in different applications. A favorable feature 
of such material is its ability to absorb large quantities of water 
from the material of high humidity. 

In this work the results of the investigation on possible 
using the starch as an absorber (hydrophilic material), are 
presented. 


EXPERIMENTAL PROCEDURE 
Material 


For the tests in question was used a sedimentary lead 
concentrate (sedimentary galena) composed mainly of PbS. 

The main component of the concentrate is the mineral 
galena whose chemical formula 1s PbS (lead sulphide). The 
sedimentary galena 1s a product of gravimetric separation of 
large mineral particles. The lead content in the sedimentary 
galena amounts to about 80%. The water content in the 
sedimentary concentrate is equal to 1-296. The sedimentary 
galena is empirically judged as a material “which may liquefy 
if shipped above the TML”. 

The following starch materials were tested (potato starch 
obtained from Potato Industry Company, Lubon, Poland): 
- cross-linked starch - diamylaceus phosphate - *Lubostat" 
- cross-linked starch - acetylic diamylaceus adipate - “AD” 
- natural starch - a granulated product 


The samples of starch materials were in a granular form. 
The experiments were conducted for samples of a concentrate 
without starch materials and of mixtures containing 99.5%, 99% 
and 98% of the concentrate and 0.5%, 1%, and 2% of starch 
materials, respectively. The used starch materials are classified 
as a product of a low environmental impact. 


Methods 


The influence of adding starch materials to the mineral 
concentrates on their parameters which determine ability for 
safe shipment by sea, was assessed on the basis of determination 
of the following parameters: 


Permeability of concentrates 

The permeability is the rate at which water under pressure 
can diffuse through the voids in the mineral concentrates. Such 
materials are permeable to water because the voids between 
the particles are inter-connected. The degree of permeability 


is characterized by the permeability coefficient k, also referred 
to as hydraulic conductivity. 

According to the classification of soils, based on their 
permeability coefficient, mineral concentrates are the materials 
with a low degree of permeability. The permeability of mineral 
concentrates depends primarily on the size and shape of grains, 
shape and arrangement of voids, void ratio, degree of saturation, 
and temperature. 

Permeability is an important parameter which determines 
the ability to liquefy the solid bulk cargoes in conditions of 
maritime transport. Liquefaction occurs first at the bottom 
layer of cargo as a result of filtration of water from the upper 
layers of cargo and for this reason its permeability should be 
determined. 


Measurement of cohesion and internal friction angle 

The estimation of cohesion and internal friction angle 
was performed in a direct shear apparatus by carrying the 
shearing with the use of the lower and upper part of displacing 
box containing the tested concentrate. In the experiment the 
samples were compacted in a dry state. The moisture content 
corresponded to the TML value estimated by means of Proctor 
C/Fagerberg Test. 

The cohesion does not occur in dry materials with pores 
entirely filled with air nor in moist materials having pores 
entirely filled with water. The cohesion increases with the 
increasing content of water and reaches its maximum values at 
a moisture approaching the TML, and then it goes down again. 
The internal friction angle reaches its minimum value when the 
moisture content is choses close to TML [13]. 

The 0.5% and 2% content of starch material in mixture 
with sedimentary galena does not significantly changes grain 
size distribution. Starch absorbed water from the mixtures at 
the amount approximately proportional to the starch material 
content in the mineral concentrate. It can be noticed that 
modified starch shows a higher solubility than natural starch. 
In general, higher TML values were observed in the case 
of testing the sedimentary galena mixed with 2% of starch 
material [14]. 

In the experiments for estimation of the permeability, 
cohesion and internal friction angle, the samples were 
compacted. The consolidation conditions in the holds were 
simulated by applying the vertical loads equal to: 0 N, 98 N, 
196 N, 294 N and 490 N, corresponding to the normal stresses 
of: 0, 1.532 - 104 N/m2, 3.0645 - 104 N/m?, 4.589 · 104 N/m2, 
7.659 - 104 N/m2, respectively. The test of samples without 
stresses corresponded to the stress conditions of cargo in 
the hold during its loading. The increasing values of normal 
stresses represented the changes in the bulk cargoes during 
their shipping by sea. 


RESULTS AND DISCUSSION 


The results of the permeability tests are presented in 
Tab. 1. 

Based on the results of the tests, the effect of different 
content of the samples on their permeability was observed. 

The maximum value of the permeability coefficient k 
was achieved for the mineral concentrate without any starch 
materials; it reached the value of 1.46 - 102 m/s. The value 
indicates that the tested concentrate may be classified as prone 
to liquefaction because its permeability coefficient significantly 
exceeds the limiting value К = 10-8 m/s [15]. 

Large spaces between particles, confirmed by the porosity 
ratio, mean that the sedimentary galena is the cargo of the 
highest permeability among the tested samples. 
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Tab.1. Results of permeability test 


Permeability coefficient [m/s] 


Sample type 


Sedimentary lead concentrate 


Sedimentary lead concentrate + 1% of starch material 
Sedimentary lead concentrate + 2% of starch material 


1.46 · 102 


Sedimentary lead concentrate + Lubostat 


1.8: 103 5.82: 104 


Sedimentary lead concentrate + AD 


2.34 · 103 1.12 · 103 


Sedimentary lead concentrate + granulated starch 


The presence of starch materials in the concentrate material 
reduces, by absorption of water, their ability to filtration. 
The ability to absorption of water is related primarily to the 
composition of starch material and the percentage of starch 
in a mixture. 

For all samples with starch material the decrease of 
permeability was found. If the content of starch material in 
the mineral concentrate increases the value of permeability 
coefficient decreases. The greatest decrease of the permeability 
coefficient was observed in the mixture containing 2% of 
Lubostat. 


Measurement of cohesion 
and internal friction angle 


The changes of cohesion in function of moisture content 
are presented in Fig. 1 through 3. 


—— sedimentary galena 


8000 — . 
——ж*— sedimentary galena + lubostat 2% 
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Fig. I. The changes of cohesion in function of moisture content for 
sedimentary galena + Lubostat 
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Fig. 2. The changes of cohesion in function of moisture content for 
sedimentary galena + AD 
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Fig. 3. The changes of cohesion in function of moisture content for 
sedimentary galena + granulated starch 


The result indicates that the starch materials significantly 
change the cohesion. The presence of the starch materials in 
mineral concentrates tends to reduce cohesion. In all tested 
samples cohesion increases with the increasing of moisture 
content and reaches its maximum with the moisture approaching 
the TML, and then it goes down. The extreme values are not 
reached at the same moisture content for all mixtures. The 
maximum values of cohesion are affected by the type of 
mixture. The result confirms presence of correlations between 
TML values of mixtures and maximum values of cohesion. 
For mixtures containing 0.5% of starch material, the greatest 
decreasing of maximum cohesion was observed. 

The changes of internal friction angle in function of 
moisture content are presented in Fig. 4 through 6. 
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Fig.4. The internal friction angle in function of moisture content for 
sedimentary galena + Lubostat 


For the sedimentary galena with starch material no 
significant changes in values of internal friction angle are 
observed. The presence of starch materials in the tested mineral 
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Fig. 5. The internal friction angle in function of moisture content for 
sedimentary galena + AD 
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Fig. 6. The internal friction angle in function of moisture content for 
sedimentary galena * granulated starch 


concentrate influences the minimum value of internal friction 
angle. In each case the values are higher than those for the 
sample without starch material. For mixtures containing 0.596 
of starch material the greatest increase of minimum internal 
friction angle was observed. The minimum values are not 
reached at the same moisture content for all samples. Content 
of starch material in mixture contributes to increase of the 
moisture content at which the maximum value of internal 
friction angle is achieved. Increasing the internal friction angle 
reduces the proneness of mineral concentrate to shifting. 


CONCLUSION 


Forall samples with starch material, a reduction ofthe degree 
of permeability was obtained. Starch materials absorbed water 
from the mixtures at the amount approximately proportional to 
the starch material content in mineral concentrates. The results 
ofthe tests indicate that the permeability influences significantly 
the properties of the concentrates, important from the point of 
view of sea transport safety. 

Water absorption and swelling the particles results in 
a reduction of space between the particles and, as the migration 
of water is thus more difficult, the concentrate capability to 
filtrate water into the lower part of cargo where the liquefaction 
process usually starts, is also reduced. 


In consequence, starch materials prevent mineral 
concentrates against drainage of water from pores between 
particles, and consequently, sliding and shifting such cargoes 
in storage. 

The comparison ofthe cohesion, internal friction angle and 
TML values confirms that for the tested mixtures the correlation 
occurs between values ofthe quantities in question. Due to the 
presence of starch materials in the mixtures the risk of passing 
mineral concentrates into the liquid state 1s lower. 
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The Ship Handling Research and Training Centre at Ilawa is owned by the Foundation for Safety of Navigation 
and Environment Protection, which is a joint venture between the Gdynia Maritime University, the Gdansk University 
of Technology and the City of Ilawa. 


Two main fields of activity of the Foundation are: 


2 Training on ship handling. Since 1980 more than 2500 ship masters and pilots from 35 countries were trained at 
Hawa Centre. The Foundation for Safety of Navigation and Environment Protection, being non-profit organisation 
Is reinvesting all spare funds in new facilities and each year to the existing facilities new models and new training 
areas were added. Existing training models each year are also modernised, that's why at present the Centre represents 
a modern facility perfectly capable to perform training on ship handling of shipmasters, pilots and tug masters. 


Research on ship's manoeuvrability. Many experimental and theoretical research programmes covering different 
problems of manoeuvrability (including human effect, harbour and waterway design) are successfully realised at 
the Centre. 


Тће Foundation possesses ISO 9001 quality certificate. 
Why training on ship handling? 


Тће safe handling of ships depends on many factors - on ship's manoeuvring characteristics, human factor (operator 
experience and skill, his behaviour in stressed situation, etc.), actual environmental conditions, and degree of water 
area restriction. 


Results of analysis of CRG (collisions, rammings and groundings) casualties show that in one third of all the 
human error is involved, and the same amount of CRG casualties is attributed to the poor controllability of ships. 
Training on ship handling is largely recommended by IMO as one of the most effective method for improving the 
safety at sea. The goal of the above training 1s to gain theoretical and practical knowledge on ship handling in a wide 
number of different situations met in practice at sea. 


For further information please contact: 
The Foundation for Safety of Navigation and Environment Protection 


Head office: Ship Handling Centre: 
36, Chrzanowskiego Street 14-200 ILAWA-KAMIONKA, POLAND 
80-278 GDANSK, POLAND tel./fax: +48 (0) 89 648 74 90 
tel./fax: +48 (0) 58 341 59 19 e-mail: office(gilawashiphandling.com.pl 

e-mail: office@portilawa.com 
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